
modeltype<-"simplified"
typename<-"stationary"
source("./auxfiles/loaddata.R")

setwd(base_dir)

setwd("./output")
obsdraws<- as.data.table(read.dta13(paste0("psid_obs",typename,".dta")))
obsdraws_cons<-as.data.table(read.dta13(paste0("psid_obs_cons",typename,".dta")))
#NEED TO CORRECT DATAMOMENTS AND INVVARCOV:
temp<-read.table(paste0("wagemomentsstationary_jump.txt"),row.names=1,fill=TRUE)
temp<-temp[rownames(temp)%in%c("prodtype_1","prodtype_2","prodtype_3","age","age2","mod","sev","old","mod_old","sev_old"),]
wageregmoments<-temp[,1]
wageregse<-as.numeric(temp[,2])
temp<-read.table(paste0("spwagemomentsstationary_jump.txt"),row.names=1,fill=TRUE)
temp<-temp[rownames(temp)%in%c("prodtype_1","prodtype_2","prodtype_3","sp_age","sp_age2","sp_mod","sp_old","sp_mod_old"),]
spwageregmoments<-temp[,1]
spwageregse<-as.numeric(temp[,2])

temp<-read.table(paste0("wagevarsstationary_jump.txt"),row.names=1,fill=TRUE)
wagevarmoments<-temp[,1]
wagevarse<-temp[,2]

# temp <-read.table(paste0("wagemeansstationary_jump.txt"),row.names=1,fill=TRUE)
# wagemeanmoments<-temp[,1]
# wagemeanse<-temp[,2]
# temp <-read.table(paste0("spwagemeansstationary_jump.txt"),row.names=1,fill=TRUE)
# spwagemeanmoments<-temp[,1]
# spwagemeanse<-temp[,2]


#invvarcov<-solve(as.matrix(read.table("~/Dropbox/DynamicDI/Low and Pistaferri/Code LP AER_Simulations/Code HL/varcovar30.inp",row.names=NULL)))
temp<-as.data.table(read.table("consmoments_all.txt"))
temp$V2<-as.numeric(as.character(temp$V2))
temp$V3<-as.numeric(as.character(temp$V3))
#dropping PT work for spouses:
temp<-temp[!grepl("sp_part_part",temp$V1),]
temp<-temp[!grepl("sp_part_part",temp$V1),]
temp<-temp[!grepl("part_part_",temp$V1)]
consmoments_all<-temp[,V2] 
consse_all<-temp[,V3] 
temp<-as.data.table(read.table("consmoments_single.txt"))
temp<-temp[!grepl("part_part_",temp$V1)]
temp$V2<-as.numeric(as.character(temp$V2))
temp$V3<-as.numeric(as.character(temp$V3))
consmoments_single<-temp[,V2]
consse_single<-temp[,V3]  
temp<-as.data.table(read.table("consmoments_married.txt"))
#dropping PT work for spouses:
temp<-temp[!grepl("sp_part_part",temp$V1),]
temp$V2<-as.numeric(as.character(temp$V2))
temp$V3<-as.numeric(as.character(temp$V3))
consmoments_married<-temp[,V2] 
consse_married<-temp[,V3] 


temp<-as.data.table(read.table("consmoments_simple_all.txt"))
temp$V2<-as.numeric(as.character(temp$V2))
temp$V3<-as.numeric(as.character(temp$V3))
#dropping PT work for spouses:
temp<-temp[!grepl("sp_part_part",temp$V1),]
temp<-temp[!grepl("part_part",temp$V1)]
consmoments_simple_all<-temp[,V2] 
consse_simple_all<-temp[,V3] 


##dropping PTwork of head and partners:
#consmoments_all<-consmoments_all[-c(4,11,14,15,18,20)]
#consse_all<-consse_all[-c(4,11,14,15,18,20)]
temp<-read.table("workmoments_agg.txt")
temp<-temp[!grepl("E=1/",rownames(temp)),]
workmoments_simple<-temp[,1]
workse_simple<-temp[,2] 


temp<-read.table("workmoments_1.txt")
temp<-temp[!grepl("E=1/",rownames(temp)),]
workmoments_mar<-temp[,1]
workse_mar<-temp[,2] 
#Dropping PT work:
#workmoments_mar<-temp[1:6,1]
#workse_mar<-temp[1:6,2] 
temp<-read.table("workmoments_0.txt")
temp<-temp[!grepl("E=1/",rownames(temp)),]
workmoments_sing<-temp[,1]
workse_sing<-temp[,2] 
# workmoments_sing<-temp[1:6,1]
# workse_sing<-temp[1:6,2] 
temp<-read.table("spworkmoments_1.txt")
#dropping PT work for spouses:
temp<-temp[!grepl("E=1/",rownames(temp)),]
spworkmoments<-temp[,1]
spworkse<-temp[,2]
#spworkmoments<-temp[1:4,1]
#spworkse<-temp[1:4,2]
temp<-read.table("stockDImoments_agg.txt")
stockDImoments_simple<-temp$EST[]
DIse_simple <- temp$SE[]

temp<-read.table("stockDImoments_work.txt")
temp<-read.table("stockDImoments.txt")
stockDImoments<-temp$EST[]
DIse <- temp$SE[]
#stockDImoments<-temp$EST[-c(1,4,7,10)]
#DIse <- temp$SE[-c(1,4,7,10)]
temp<-read.table("compDImoments_agg.txt")
compDImoments_simple<-temp$EST[]
compDIse_simple<-temp$SE[]

temp<-read.table("compDImoments.txt")
compDImoments<-temp$EST[]
compDIse<-temp$SE[]
#compDImoments<-temp$EST[-c(5,6,11,12)]
#compDIse<-temp$SE[-c(5,6,11,12)]

#setwd("../eventstudy/PSID_eventstudy")
rform_events_mod<-as.data.table(read.csv("../eventstudy/PSID_eventstudy/mod_eventmeans_table.csv"))
rform_events_sev<-as.data.table(read.csv("../eventstudy/PSID_eventstudy/sev_eventmeans_table.csv"))
eventDImoments_single<-c(unlist(rform_events_mod[X%in%c("as.factor(married)0"),'Estimate']),
                         unlist(rform_events_sev[X%in%c("as.factor(married)0"),'Estimate']))
eventDIse_single<-c(unlist(rform_events_mod[X%in%c("as.factor(married)0"),'Std..Error']),
                    unlist(rform_events_sev[X%in%c("as.factor(married)0"),'Std..Error']))
eventDImoments_married<-c(unlist(rform_events_mod[X%in%c("as.factor(married)1"),'Estimate']),
                          unlist(rform_events_sev[X%in%c("as.factor(married)1"),'Estimate']))
eventDIse_married<-c(unlist(rform_events_mod[X%in%c("as.factor(married)1"),'Std..Error']),
                     unlist(rform_events_sev[X%in%c("as.factor(married)1"),'Std..Error']))

eventDImoments<-c(eventDImoments_single,eventDImoments_married)
eventDIse<-c(eventDIse_single,eventDIse_married)

rform_events_mod<-as.data.table(read.csv("../eventstudy/PSID_eventstudy/mod_agg_eventmeans_table.csv"))
rform_events_sev<-as.data.table(read.csv("../eventstudy/PSID_eventstudy/sev_agg_eventmeans_table.csv"))
eventDImoments_simple<-c(unlist(rform_events_mod[X%in%c("Estimate"),'x']),
                         unlist(rform_events_sev[X%in%c("Estimate"),'x']))
eventDIse_simple<-c(unlist(rform_events_mod[X%in%c("Std. Error"),'x']),
                    unlist(rform_events_sev[X%in%c("Std. Error"),'x']))

#setwd("../../output")

# temp1<-read.table("eventheadmodDImoments_ddbal.txt")
# temp2<-read.table("eventheadsevDImoments_ddbal.txt")
# temp1[,"V3"]<-as.character(temp1[,"V3"])
# temp1[,"V3"]<-as.numeric(gsub("\\(","",gsub("\\)","",temp1[,"V3"])))
# temp2[,"V3"]<-as.character(temp2[,"V3"])
# temp2[,"V3"]<-as.numeric(gsub("\\(","",gsub("\\)","",temp2[,"V3"])))
# eventDImoments<-c(temp1[1,'V2'],temp2[1,'V2'],temp1[2,'V2'],temp2[2,'V2'])
# eventDIse<-c(temp1[1,'V3'],temp2[1,'V3'],temp1[2,'V3'],temp2[2,'V3'])

temp<-read.table("flowDImoments_agg.txt")
flowDImoments_simple<-temp$EST[]
flowDIse_simple<- temp$SE[]


temp<-read.table("flowDImomentsbymar.txt")
flowDImoments<-temp$EST[]
flowDIse<- temp$SE[]
#flowDImoments<-temp$EST[-c(1,4)]
#flowDIse<- temp$SE[-c(1,4)]

temp<-read.table("contflowDImoments_agg.txt")
contflowDImoments_simple<-temp$EST
names(contflowDImoments_simple)<-rownames(temp)
contflowDIse_simple<- temp$SE
names(contflowDIse_simple)<-rownames(temp)

temp<-read.table("contflowDImomentsbymar.txt")
contflowDImoments<-temp$EST
names(contflowDImoments)<-rownames(temp)
contflowDIse<- temp$SE
names(contflowDIse)<-rownames(temp)


#dropping spousal moments:
#consmoments_all<-consmoments_all[-c(12,13,14)]
#consse_all<-consse_all[-c(12,13,14)]
#spworkmoments<-NULL
#spworkse<-NULL

temp<-read.table("jointworkmoments.txt")
jointworkmoments<-temp$EST #(no work, no spwork), (no work, sp work), (work, no spwork), (work, spwork)

temp<-read.table("jointworkpartmoments.txt")
jointworkpartmoments<-temp$EST #(no work, no spwork), (no work, sp work), (work, no spwork), (work, spwork)
#Including event moments instead of flow and composition moments:
datamoments<-c(wageregmoments,spwageregmoments,wagevarmoments,
               consmoments_simple_all,workmoments_simple,spworkmoments, 
               #stockDImoments,
               compDImoments_simple,
               flowDImoments_simple,
               #               contflowDImoments
               eventDImoments_simple
)

datases<-c(wageregse,spwageregse,wagevarse,
               consse_simple_all,workse_simple,spworkse, 
               #stockDImoments,
               compDIse_simple,
               flowDIse_simple,
               #               contflowDImoments
               eventDIse_simple
)

wagevarnames<-c("Variance","Autocovariance","Spouse Variance","Spouse Autocovariance","Intra-Household Covariance")
eventDInames <- c("DI after Moderate Shock","DI after Severe Shock")
plotwagenames<-c("Low Type","Mid Type","High Type","Sev. Dis.","Mod. Dis.","Age","Age Sq.","Age \u2265 45","Age \u2265 45 and Sev. Dis.","Age \u2265 45 and Mod. Dis.")
plotspwagenames<-c("Low Type","Mid Type","High Type","Disabled","Age","Age Sq.","Age \u2265 45","Age \u2265 45 and Disabled")
plotconsnames_all<-c("Disability","Work","DI Receipt","Work X Disability","Spousal Disability","Spousal Work","Spousal Work X Disability")
plotworknames<-c("Healthy, Age < 45","Healthy, Age \u2265 45","Mod. Dis., Age < 45", "Mod. Dis., Age \u2265 45","Sev. Dis., Age < 45","Sev. Dis., Age \u2265 45")
plotspworknames <- c("Healthy, Age < 45","Healthy, Age \u2265 45", "Disabled, Age < 45", "Disabled, Age \u2265 45")
plotcompDInames <- c("Sev. Dis., Age < 45","Sev. Dis., Age \u2265 45", "Mod. Dis., Age < 45","Mod. Dis., Age \u2265 45","Healthy, Age < 45","Healthy, Age \u2265 45")
plotflowDInames <- c("Healthy, Age < 45","Healthy, Age \u2265 45","Mod. Dis., Age < 45","Mod. Dis., Age \u2265 45","Sev. Dis., Age < 45","Sev. Dis., Age \u2265 45")


setwd(base_dir)

path<-"./modeloutput/"

#Matrices for Counterfactual SEs:
#First three rows are gradients for wtp, fiscal cost, mvpf, employment. 11th row is DI Rolls:
# Jmat<-read.csv(paste0("./modeloutput/newguess_jacobian_",modeltype,".csv"),header=FALSE)
# Cfacts<-list(
# SSgenerosity=read.csv(paste0("./modeloutput/gradient_ss_",modeltype,".csv"),header=FALSE),
# FSgenerosity=read.csv(paste0("./modeloutput/gradient_fs_",modeltype,".csv"),header=FALSE),
# stringency=read.csv(paste0("./modeloutput/gradient_stringency_",modeltype,".csv"),header=FALSE),
# singstringency=read.csv(paste0("./modeloutput/gradient_singstringency_",modeltype,".csv"),header=FALSE),
# marstringency=read.csv(paste0("./modeloutput/gradient_marstringency_",modeltype,".csv"),header=FALSE),
# EarlyDI=read.csv(paste0("./modeloutput/gradient_earlyDI_",modeltype,".csv"),header=FALSE),
# MeansTest_1=read.csv(paste0("./modeloutput/gradient_meanstest_",modeltype,".csv"),header=FALSE),
# SSgenerositymeanstest_1=read.csv(paste0("./modeloutput/gradient_ss_meanstest_",modeltype,".csv"),header=FALSE),
# FSgenerositymeanstest_1=read.csv(paste0("./modeloutput/gradient_fs_meanstest_",modeltype,".csv"),header=FALSE),
# stringencymeanstest_1=read.csv(paste0("./modeloutput/gradient_stringency_meanstest_",modeltype,".csv"),header=FALSE),
# EarlyDImeanstest_1=read.csv(paste0("./modeloutput/gradient_earlydi_meanstest_",modeltype,".csv"),header=FALSE)
# )
# vcovnames<-c("wtp","constax","mvpf",
#              "work","work.m.1","work.m.2","work.m.3","work.s.1","work.s.2","work.s.3",
#              "DI","DI.m.1","DI.m.2","DI.m.3","DI.s.1","DI.s.2","DI.s.3")
# for(v in 1:length(Cfacts)){
#   rownames(Cfacts[[v]])<-vcovnames
# }

paramvars<-tryCatch(read.csv(paste0("./modeloutput/newguess_varcov_",modeltype,".csv"),header=FALSE),error=function(x) matrix(NA,length(newguess[which(upper!=lower)]),length(newguess[which(upper!=lower)])))
Lambda<-tryCatch(read.csv(paste0("./modeloutput/newguess_sensitivity_",modeltype,".csv"),header=FALSE),error=function(x) matrix(NA,length(newguess[which(upper!=lower)]),length(newguess[which(upper!=lower)])))
paramses<-diag(as.matrix(paramvars))
paramses<-paramses^0.5

temp<-tryCatch(unique(as.data.table(read.csv(file=paste0("./modeloutput/bootstrap/bootstrap_","CV","_",modeltype,".csv"),header=FALSE))),error=function(x) NA)
temp<-temp[1:200,] #just so I can say we used exactly 200 replications.
temp<-temp[,1:8]
names(temp)<-c("theta_sev","sptheta","eta","speta","eta_sev","speta_mod","F","spF")
cv_vars<-temp
cv_vars<-apply(cv_vars,MARGIN=2,FUN=function(x)as.numeric(gsub("\\\\\\$","",x)))
rm(temp)

names(paramses)<-names(newguess[which(upper!=lower)])
rownames(Lambda)<-names(newguess[which(upper!=lower)])
momentnames<-c(plotwagenames,plotspwagenames,wagevarnames,plotconsnames_all,plotworknames,plotspworknames,plotcompDInames,plotflowDInames,eventDInames)
momentgroups<-c(rep("wage",length(plotwagenames)),
                rep("spwage",length(plotspwagenames)),
                rep("wagevar",length(wagevarnames)),
                rep("cons",length(plotconsnames_all)),
                rep("work",length(plotworknames)),
                rep("spwork",length(plotspworknames)),
                rep("DI",length(plotcompDInames)),
                rep("DI",length(plotflowDInames)),
                rep("DI",length(eventDInames)))

colnames(Lambda)<-momentnames
Lambda[c("theta_sev","eta_sev","eta"),c(#4,5,9,10, #realized wages
                                              24:27, #consumption
                                              31:36)] #work

Lambda_scaled<-sapply(1:length(datases),function(x) Lambda[,x])
colnames(Lambda_scaled)<-colnames(Lambda)
rownames(Lambda_scaled)<-rownames(Lambda)
#This just favors most precisely-estimated moments...
#excluding healthy DI flows -- these are zero.
tab_sens<-sapply(1:nrow(Lambda_scaled[,]),FUN=function(x){
  y<-data.table(val=abs(unlist(Lambda_scaled[x,])),name=momentgroups)
  y<-y[,max(val),by=name]
  z<-y$V1/y$V1[1]
# if(x %in% 1:7)  z<-y$V1/y$V1[1]
# if(x %in% 8:13)  z<-y$V1/y$V1[2]
# if(x %in% 14:18)  z<-y$V1/y$V1[3]
# if(x %in% 19:24)  z<-y$V1/y$V1[4]
# if(x %in% c(25,27))  z<-y$V1/y$V1[5]
# if(x %in% c(26,28))  z<-y$V1/y$V1[6]
# if(x %in% 29:34)  z<-y$V1/y$V1[8]
  names(z)<-y$name
  
  return(z)
} )
colnames(tab_sens)<-rownames(Lambda)
tab_sens<-t(tab_sens)
tab_sens<-tab_sens[,-c(2,3,6)]

printparams<-TR(c("","Wage","Consumption","Work","DI"))+
  midrulep(list(c(2,5))) +
  TR("Panel A: Preference Parameters",cspan=5,surround = "\\textit{%s}")+
  TR(c("$\\theta^1$"))%:%TR(c(tab_sens[rownames(tab_sens)=='theta_sev',]),cspan=c(1),dec=3)+
  TR(c("$\\eta^1$"))%:%TR(c(tab_sens[rownames(tab_sens)=='eta',]),cspan=c(1),dec=3)+
  TR(c("$\\eta^{1,dis}$"))%:%TR(c(tab_sens[rownames(tab_sens)=='eta_sev',]),cspan=c(1),dec=3)+
#  TR(c("$\\theta^{2}$"))%:%TR(c(tab_sens[rownames(tab_sens)=='sptheta',]),cspan=c(1),dec=3)+
#  TR(c("$\\eta^{2}$"))%:%TR(c(tab_sens[rownames(tab_sens)=='speta',]),cspan=c(1),dec=3)+
#  TR(c("$\\eta^{2,dis}$"))%:%TR(c(tab_sens[rownames(tab_sens)=='speta_mod',]),cspan=c(1),dec=3)+
  vspace(5)+
  TR("Panel B: Frictions and Fixed Costs",cspan=5,surround = "\\textit{%s}")+
  TR(c("$\\delta^{1}$"))%:%TR(c(tab_sens[rownames(tab_sens)=='delta',]),cspan=c(1),dec=3)+
#  TR(c("$\\delta^{2}$"))%:%TR(c(tab_sens[rownames(tab_sens)=='spdelta',]),cspan=c(1),dec=3)+
  TR(c("$F^1$"))%:%TR(c(tab_sens[rownames(tab_sens)=='F3',]),cspan=c(1),dec=3)+
#  TR(c("$F^2$"))%:%TR(c(tab_sens[rownames(tab_sens)=='spF2',]),cspan=c(1),dec=3)+
  vspace(5)+
  TR("Panel C: Allowance Probabilities",cspan=5,surround = "\\textit{%s}")+
  TR(c("$\\pi_0^{young}$"))%:%TR(c(tab_sens[rownames(tab_sens)=='pi0.young.married',]),cspan=c(1),dec=3)+
  TR(c("$\\pi_1^{young}$"))%:%TR(c(tab_sens[rownames(tab_sens)=='pi1.young.married',]),cspan=c(1),dec=3)+
  TR(c("$\\pi_2^{young}$"))%:%TR(c(tab_sens[rownames(tab_sens)=='pi2.young.married',]),cspan=c(1),dec=3)+
  TR(c("$\\pi_0^{old}$"))%:%TR(c(tab_sens[rownames(tab_sens)=='pi0.old.married',]),cspan=c(1),dec=3)+
  TR(c("$\\pi_1^{old}$"))%:%TR(c(tab_sens[rownames(tab_sens)=='pi1.old.married',]),cspan=c(1),dec=3)+
  TR(c("$\\pi_2^{old}$"))%:%TR(c(tab_sens[rownames(tab_sens)=='pi2.old.married',]),cspan=c(1),dec=3)+
  TR("Panel D: Wage Coefficients",cspan=c(5), surround = "\\textit{%s}")+
  TR(c("$Age$"))%:%TR(c(tab_sens[rownames(tab_sens)=='age',]),cspan=c(1),dec=3)+
  TR(c("$Age^2$"))%:%TR(c(tab_sens[rownames(tab_sens)=='age2',]),cspan=c(1),dec=3)+
  TR(c("$H^{mod}$"))%:%TR(c(tab_sens[rownames(tab_sens)=='mod',]),cspan=c(1),dec=3)+
  TR(c("$H^{sev}$"))%:%TR(c(tab_sens[rownames(tab_sens)=='sev',]),cspan=c(1),dec=3)+
  TR(c("$f_1$"))%:%TR(c(tab_sens[rownames(tab_sens)=='Type1',]),cspan=c(1),dec=3)+
  TR(c("$f_2$"))%:%TR(c(tab_sens[rownames(tab_sens)=='Type2',]),cspan=c(1),dec=3)+
  TR(c("$f_3$"))%:%TR(c(tab_sens[rownames(tab_sens)=='Type3',]),cspan=c(1),dec=3)
#  vspace(5)+
  # TR("Panel E: Spousal Wage Coefficients",cspan=c(5), surround = "\\textit{%s}")+
  # TR(c("$Age$"))%:%TR(c(tab_sens[rownames(tab_sens)=='spage',]),cspan=c(1),dec=3)+
  # TR(c("$Age^2$"))%:%TR(c(tab_sens[rownames(tab_sens)=='spage2',]),cspan=c(1),dec=3)+
  # TR(c("$H^{mod}$"))%:%TR(c(tab_sens[rownames(tab_sens)=='spmod',]),cspan=c(1),dec=3)+
  # TR(c("$f_1$"))%:%TR(c(tab_sens[rownames(tab_sens)=='spType1',]),cspan=c(1),dec=3)+
  # TR(c("$f_2$"))%:%TR(c(tab_sens[rownames(tab_sens)=='spType2',]),cspan=c(1),dec=3)+
  # TR(c("$f_3$"))%:%TR(c(tab_sens[rownames(tab_sens)=='spType3',]),cspan=c(1),dec=3)



TS(printparams,file=paste0("paramsensitivities_all_",modeltype),
   output_path=path,pretty_rules=T,
   header=c("r",rep("c",4)))

#source("../StructuralModel/Welfare_ParallelSpouseC.R")
#setThreadOptions(numThreads=20,stackSize=1e07)

#sourceCpp("~/Dropbox/DynamicDI/StructuralModel/cpp/Valsim_corrwage_wpolicy.cpp", rebuild=TRUE)

#source("../StructuralModel/auxfiles/loaddata.R")

if(health_simpleprefs == TRUE) {
  Lambda[c("theta_sev","eta_sev","F3"),] <- Lambda[c("theta_sev","eta_sev","F3"),]/2
}


plot_lambda(Lambda, path=path, suffix=modeltype, nocomp=nocomp,noflow=noflow,noeventDI=noeventDI,
            nocontflow=nocontflow, nostockDI=nostockDI)

reformnames<-c("EarlyDI","stringency","singstringency","SSgenerosity","FSgenerosity")
reformbalnames<-c("Temp. DI",  "DI Stringency (%)",  "DI Stringency, single only (%)", "DI Generosity (%)", "Welfare Generosity (%)")
reformnames_meanstest<-c("EarlyDImeanstest_1","stringencymeanstest_1","singstringencymeanstest_1","SSgenerositymeanstest_1","FSgenerositymeanstest_1")
reformnames_meanstestonly<-c("MeansTest_1","MeansTest_20000","MeansTest_40000","MeansTest_90000")
reformnames_meanstestnoSSDI<-c("EarlyDImeanstestnoSSDI","stringencymeanstestnoSSDI","singstringencymeanstestnoSSDI","SSgenerositymeanstestnoSSDI","FSgenerositymeanstestnoSSDI")
#WorkDI WorkApply
reformlabels<-c("Temp. DI",
                "DI Leniency","DI Leniency (single only)","DI Generosity",
                "Food Stamp \n Generosity")
reformlabels_meanstest<-c("Temp. DI",
                "DI Leniency","DI Leniency (single only)","DI Generosity", "Food Stamp \n Generosity")
reformlabels_meanstestonly<-c("No Savings",
                          "Savings $\\le 20,000$","Savings $\\le 40,000$", "Savings $\\le 90,000$")
reformlabels_meanstestnoSSDI<-c("Temp. DI",
                          "DI Leniency","DI Leniency (single only)","DI Generosity", "Food Stamp \n Generosity")


simdata<-readRDS("./modeloutput/simdata.RDS")
simdata<-cleansim(simdata,select=select,selectC=selectC)

uflow<-((simdata[age<62,mean(flow)]*(1-params$gamma))^(1/(1-params$gamma)))

#"Allow  App. Work",
if(modeltype=="simplified"){
printparams<-TR(c("","Male Worker","Female Spouse"),cspan=c(1,4,4))+
  TR(c("","Point Estimate","Cost","Point Estimate","Cost"),cspan=c(1,2,2,2,2))+
  TR(c("","Coef","St. Error","Avg. C.V.","St. Error","Coef","St. Error","Avg. C.V.","St. Error"))+
  midrulep(list(c(2,3),c(4,5),c(6,7),c(8,9))) +
  TR("Panel A: Preferences",cspan=6,surround = "\\textit{%s}")+
  TR(c("$\\theta$"))%:%TR(c(newguess[c("theta_sev")]/2),cspan=c(1),dec=3)%:%TR(c(paramses[c("theta_sev")]/2),cspan=c(1),se=TRUE,dec=3)%:%
  TR(cost(newguess[c("theta_sev")]/2,uflow))%:% TR(paste0("(",as.character(round(var(cv_vars[,"theta_sev"])^0.5)),")"))%:%
  TR(newguess["sptheta"],cspan=1,dec=3)%:%TR(paramses["sptheta"],cspan=1,se=TRUE,dec=3)%:%
  TR(cost(newguess[c("sptheta")],uflow))%:%TR(paste0("(",as.character(round(var(cv_vars[,"sptheta"])^0.5)),")"))+
  TR(c("$\\eta$"))%:%TR(c(newguess["eta"]),cspan=c(1),dec=3)%:%TR(paramses["eta"],cspan=c(1),se=TRUE,dec=3)%:%
  TR(cost(newguess[c("eta")],uflow))%:% TR(paste0("(",as.character(round(var(cv_vars[,"eta"])^0.5)),")"))%:%
  TR(newguess["speta"],cspan=1,dec=3)%:%TR(paramses["speta"],cspan=1,se=TRUE,dec=3)%:%
  TR(cost(newguess[c("speta")],uflow))%:% TR(paste0("(",as.character(round(var(cv_vars[,"speta"])^0.5)),")"))+
  TR("$\\eta^{dis}$")%:%TR(newguess[c("eta_sev")]/2,cspan=c(1),dec=3)%:%TR(c(paramses[c("eta_sev")]/2),cspan=c(1),dec=3,se=TRUE)%:%
  TR(cost(newguess[c("eta_sev")]/2,uflow))%:% TR(paste0("(",as.character(round(var(cv_vars[,"eta_sev"])^0.5)),")"))%:%
  TR(c(newguess["speta_mod"]),cspan=c(1),dec=3)%:%TR(paramses["speta_mod"],cspan=c(1),se=TRUE,dec=3)%:%
  TR(cost(newguess[c("speta_mod")],uflow))%:% TR(paste0("(",as.character(round(var(cv_vars[,"speta_mod"])^0.5)),")"))+
  vspace(5)+
  TR("Panel B: Frictions and Fixed Costs",cspan=6,surround = "\\textit{%s}")+
  TR(c("$\\delta$"))%:%TR(c(newguess["delta"]),cspan=c(1),dec=3)%:%TR(paramses["delta"],cspan=c(1),se=TRUE,dec=3)%:%TR(c(""))%:%TR(c(""))%:%
  TR(newguess["spdelta"],cspan=1,dec=3)%:%TR(paramses["spdelta"],cspan=1,se=TRUE,dec=3)%:%TR(c(""))+
  TR(c("$F$"))%:%TR(c(newguess["F3"]/2),cspan=c(1),dec=3)%:%TR(paramses["F3"]/2,cspan=c(1),se=TRUE,dec=3)%:%
  TR(paste0("\\$",as.character(round(param_input["F3"]/2))))%:%
  TR(paste0("(",as.character(round(var(cv_vars[,"F"])^0.5)),")"))%:%
  TR(newguess["spF2"],cspan=1,dec=3)%:%TR(paramses["spF2"],cspan=1,se=TRUE,dec=3)%:%
  TR(paste0("\\$",as.character(round(param_input["spF2"]))))%:%
  TR(paste0("(",as.character(round(var(cv_vars[,"spF"])^0.5)),")"))+
  vspace(5)+
  TR("Panel C: Allowance Probabilities",cspan=6,surround = "\\textit{%s}")+
  TR(c("$\\pi_0^{young}$"))%:%TR(c(newguess["pi0.young.married"]),cspan=c(1),dec=3)%:%TR(paramses["pi0.young.married"],cspan=c(1),se=TRUE,dec=3)%:%TR(c(""))%:%TR(c(""))%:%
  TR("---",cspan=2)%:%TR(c(""))%:%TR(c(""))+
  TR(c("$\\pi_1^{young}$"))%:%TR(c(newguess["pi1.young.married"]),cspan=c(1),dec=3)%:%TR(paramses["pi1.young.married"],cspan=c(1),se=TRUE,dec=3)%:%TR(c(""))%:%TR(c(""))%:%
  TR("---",cspan=2)%:%TR(c(""))%:%TR(c(""))+
  TR(c("$\\pi_2^{young}$"))%:%TR(c(newguess["pi2.young.married"]),cspan=c(1),dec=3)%:%TR(paramses["pi2.young.married"],cspan=c(1),se=TRUE,dec=3)%:%TR(c(""))%:%TR(c(""))%:%
  TR("---",cspan=2)%:%TR(c(""))%:%TR(c(""))+
  TR(c("$\\pi_0^{old}$"))%:%TR(c(newguess["pi0.old.married"]),cspan=c(1),dec=3)%:%TR(paramses["pi0.old.married"],cspan=c(1),se=TRUE,dec=3)%:%TR(c(""))%:%TR(c(""))%:%
  TR("---",cspan=2)%:%TR(c(""))%:%TR(c(""))+
  TR(c("$\\pi_1^{old}$"))%:%TR(c(newguess["pi1.old.married"]),cspan=c(1),dec=3)%:%TR(paramses["pi1.old.married"],cspan=c(1),se=TRUE,dec=3)%:%TR(c(""))%:%TR(c(""))%:%
  TR("---",cspan=2)%:%TR(c(""))%:%TR(c(""))+
  TR(c("$\\pi_2^{old}$"))%:%TR(c(newguess["pi2.old.married"]),cspan=c(1),dec=3)%:%TR(paramses["pi2.old.married"],cspan=c(1),se=TRUE,dec=3)%:%TR(c(""))%:%TR(c(""))%:%
  TR("---",cspan=2)%:%TR(c(""))%:%TR(c(""))+
  TR("Panel D: Wage Coefficients",cspan=c(6), surround = "\\textit{%s}")+
  TR(c("Age"))%:%TR(c(newguess[c("age")]),cspan=c(1),dec=3)%:%TR(c(paramses[c("age")]),cspan=c(1),se=TRUE,dec=3)%:%TR(c(""))%:%TR(c(""))%:%
  TR(newguess["spage"],cspan=1,dec=3)%:%TR(paramses["spage"],cspan=1,se=TRUE,dec=3)%:%TR(c(""))%:%TR(c(""))+
  TR(c("Age$^2$"))%:%TR(c(newguess["age2"]),cspan=c(1),dec=3)%:%TR(paramses["age2"],cspan=c(1),se=TRUE,dec=3)%:%TR(c(""))%:%TR(c(""))%:%
  TR(newguess["spage2"],cspan=1,dec=3)%:%TR(paramses["spage2"],cspan=1,se=TRUE,dec=3)%:%TR(c(""))%:%TR(c(""))+
  TR("$H^{mod}$")%:%TR(newguess[c("mod")],cspan=c(1),dec=3)%:%TR(c(paramses[c("mod")]),cspan=c(1),dec=3,se=TRUE)%:%TR(c(""))%:%TR(c(""))%:%
  TR(c(newguess["spmod"]),cspan=c(1),dec=3)%:%TR(paramses["spmod"],cspan=c(1),se=TRUE,dec=3)%:%TR(c(""))%:%TR(c(""))+
  TR("$H^{sev}$")%:%TR(newguess[c("sev")],cspan=c(1),dec=3)%:%TR(c(paramses[c("sev")]),cspan=c(1),dec=3,se=TRUE)%:%TR(c(""))%:%TR(c(""))%:%
  TR("---",cspan=2)%:%TR(c(""))%:%TR(c(""))+
  TR("$f_1$")%:%TR(newguess[c("Type1")],cspan=c(1),dec=3)%:%TR(c(paramses[c("Type1")]),cspan=c(1),dec=3,se=TRUE)%:%TR(c(""))%:%TR(c(""))%:%
  TR(c(newguess["spType1"]),cspan=c(1),dec=3)%:%TR(paramses["spType1"],cspan=c(1),se=TRUE,dec=3)%:%TR(c(""))%:%TR(c(""))+
  TR("$f_2$")%:%TR(newguess[c("Type2")],cspan=c(1),dec=3)%:%TR(c(paramses[c("Type2")]),cspan=c(1),dec=3,se=TRUE)%:%TR(c(""))%:%TR(c(""))%:%
  TR(c(newguess["spType2"]),cspan=c(1),dec=3)%:%TR(paramses["spType2"],cspan=c(1),se=TRUE,dec=3)%:%TR(c(""))%:%TR(c(""))+
  TR("$f_3$")%:%TR(newguess[c("Type3")],cspan=c(1),dec=3)%:%TR(c(paramses[c("Type3")]),cspan=c(1),dec=3,se=TRUE)%:%TR(c(""))%:%TR(c(""))%:%
  TR(c(newguess["spType3"]),cspan=c(1),dec=3)%:%TR(paramses["spType3"],cspan=c(1),se=TRUE,dec=3)%:%TR(c(""))%:%TR(c(""))+
  vspace(5)+
  TR("Panel E: Wage Variance Parameters",cspan=c(6), surround = "\\textit{%s}")+
  TR("$\\sigma_\\epsilon$")%:%TR(c(newguess["noisevar"]^0.5),dec=4)%:%TR(c(paramses["noisevar"]^0.5), dec=4,se=TRUE)%:%TR(c(""))%:%TR(c(""))%:%
  TR(c(newguess["spnoisevar"]^0.5),dec=4)%:%TR(c(paramses["spnoisevar"]^0.5), dec=4,se=TRUE)%:%TR(c(""))%:%TR(c(""))+
  TR("$\\sigma_\\xi$")%:%TR(c(newguess["wagevar"]^0.5),dec=4)%:%TR(c(paramses["wagevar"]^0.5), dec=4,se=TRUE)%:%TR(c(""))%:%TR(c(""))%:%
  TR(c(newguess["spwagevar"]^0.5),dec=4)%:%TR(c(paramses["spwagevar"]^0.5), dec=4,se=TRUE)%:%TR(c(""))%:%TR(c(""))+
  TR("$\\rho_\\xi$")%:%TR(c(newguess["wagecorr"]),dec=4)%:%TR(c(paramses["wagecorr"]),dec=4,se=TRUE)%:%TR(c(""))%:%TR(c(""))



TS(printparams,file=paste0("paramvals_all_",modeltype),
   output_path=path,pretty_rules=T,
   header=c("r",rep("c",8)))
} else{
  printparams<-TR(c("","Single","Married","Spouse"),cspan=c(1,2,2,2))+
    TR(c("","Coef","St. Error","Coef","St. Error","Coef","St. Error"))+
    midrulep(list(c(2,3),c(4,5),c(6,7))) +
    TR("Panel A: Preferences",cspan=7,surround = "\\textit{%s}")+
    TR(c("$\\mu$"))%:%TR(c("---"),cspan=c(2))%:%TR(c("---"),cspan=c(2))%:%
    TR(c(newguess[c("spmu")]),cspan=c(1),dec=3)%:%TR(paramses[c("spmu")],se=TRUE,cspan=c(1),dec=3)+
    TR(c("$\\theta^{mod}$"))%:%TR(c(newguess["theta_mod"]),cspan=c(1),dec=3)%:%TR(paramses["theta_mod"],cspan=c(1),se=TRUE,dec=3)%:%
    TR("(same as single)",cspan=2)%:%TR(newguess["sptheta"],cspan=1,dec=3)%:%TR(paramses["sptheta"],cspan=1,se=TRUE,dec=3)+
    TR(c("$\\theta^{sev}$"))%:%TR(c(newguess[c("theta_sev")]),cspan=c(1),dec=3)%:%TR(c(paramses[c("theta_sev")]),cspan=c(1),se=TRUE)%:%
    TR("(same as single)",cspan=2)%:%TR(c("---"),cspan=c(2))+
    TR(c("$\\eta$"))%:%TR(newguess["eta_single"],cspan=c(1),dec=3)%:%TR(paramses["eta_single"],cspan=c(1),se=TRUE,dec=3)%:%
    TR(c(newguess["eta"]),cspan=c(1),dec=3)%:%TR(paramses["eta"],cspan=c(1),se=TRUE,dec=3)%:%
    TR(newguess["speta"],cspan=1,dec=3)%:%TR(paramses["sptheta"],cspan=1,se=TRUE,dec=3)+
    TR(c("$\\eta^{mod}$"))%:%TR(newguess[c("eta_mod")],cspan=c(1),dec=3)%:%TR(c(paramses[c("eta_mod")]),cspan=c(1),se=TRUE,dec=3)%:%
    TR("(same as single)",cspan=2)%:%TR(c(newguess["speta_mod"]),cspan=c(1),dec=3)%:%TR(paramses["speta_mod"],cspan=c(1),se=TRUE,dec=3)+
    TR("$\\eta^{sev}$")%:%TR(newguess[c("eta_sev")],cspan=c(1),dec=3)%:%TR(c(paramses[c("eta_sev")]),cspan=c(1),dec=3,se=TRUE)%:%
    TR("(same as single)",cspan=2)%:%TR(c("---"),cspan=c(2))+
    vspace(5)+
    TR("Panel B: Frictions and Fixed Costs",cspan=7,surround = "\\textit{%s}")+
    TR(c("$\\delta$"))%:%TR(newguess["delta_single"],cspan=c(1),dec=3)%:%TR(paramses["delta_single"],cspan=c(1),se=TRUE,dec=3)%:%
    TR(c(newguess["delta"]),cspan=c(1),dec=3)%:%TR(paramses["delta"],cspan=c(1),se=TRUE,dec=3)%:%
    TR(newguess["spdelta"],cspan=1,dec=3)%:%TR(paramses["spdelta"],cspan=1,se=TRUE,dec=3)+
    TR(c("$F^{mod}$"))%:%TR(newguess["F2_single"],cspan=c(1),dec=3)%:%TR(paramses["F2_single"],cspan=c(1),se=TRUE,dec=3)%:%
    TR(c(newguess["F2"]),cspan=c(1),dec=3)%:%TR(paramses["F2"],cspan=c(1),se=TRUE,dec=3)%:%
    TR(newguess["spF2"],cspan=1,dec=3)%:%TR(paramses["spF2"],cspan=1,se=TRUE,dec=3)+
    TR(c("$F^{sev}$"))%:%TR(newguess["F3_single"],cspan=c(1),dec=3)%:%TR(paramses["F3_single"],cspan=c(1),se=TRUE,dec=3)%:%
    TR(c(newguess["F3"]),cspan=c(1),dec=3)%:%TR(paramses["F3"],cspan=c(1),se=TRUE,dec=3)%:%
    TR("---",cspan=2)+
    TR("$F^{old}$")%:%TR(newguess["Fold_single"],cspan=c(1),dec=3)%:%TR(paramses["Fold_single"],cspan=c(1),se=TRUE,dec=3)%:%
    TR(c(newguess["Fold"]),cspan=c(1),dec=3)%:%TR(paramses["Fold"],cspan=c(1),se=TRUE,dec=3)%:%
    TR(newguess["spFold"],cspan=1,dec=3)%:%TR(paramses["spFold"],cspan=1,se=TRUE,dec=3)+
    vspace(5)+
    TR("Panel C: Allowance Probabilities",cspan=7,surround = "\\textit{%s}")+
    TR(c("$\\pi_0^{young}$"))%:%TR(newguess["pi0.young.single"],cspan=c(1),dec=3)%:%TR(paramses["pi0.young.single"],cspan=c(1),se=TRUE,dec=3)%:%
    TR(c(newguess["pi0.young.married"]),cspan=c(1),dec=3)%:%TR(paramses["pi0.young.married"],cspan=c(1),se=TRUE,dec=3)%:%
    TR("---",cspan=2)+
    TR(c("$\\pi_1^{young}$"))%:%TR(newguess["pi1.young.single"],cspan=c(1),dec=3)%:%TR(paramses["pi1.young.single"],cspan=c(1),se=TRUE,dec=3)%:%
    TR(c(newguess["pi1.young.married"]),cspan=c(1),dec=3)%:%TR(paramses["pi1.young.married"],cspan=c(1),se=TRUE,dec=3)%:%
    TR("---",cspan=2)+
    TR(c("$\\pi_2^{young}$"))%:%TR(newguess["pi2.young.single"],cspan=c(1),dec=3)%:%TR(paramses["pi2.young.single"],cspan=c(1),se=TRUE,dec=3)%:%
    TR(c(newguess["pi2.young.married"]),cspan=c(1),dec=3)%:%TR(paramses["pi2.young.married"],cspan=c(1),se=TRUE,dec=3)%:%
    TR("---",cspan=2)+
    TR(c("$\\pi_0^{old}$"))%:%TR(newguess["pi0.old.single"],cspan=c(1),dec=3)%:%TR(paramses["pi0.old.single"],cspan=c(1),se=TRUE,dec=3)%:%
    TR(c(newguess["pi0.old.married"]),cspan=c(1),dec=3)%:%TR(paramses["pi0.old.married"],cspan=c(1),se=TRUE,dec=3)%:%
    TR("---",cspan=2)+
    TR(c("$\\pi_1^{old}$"))%:%TR(newguess["pi1.old.single"],cspan=c(1),dec=3)%:%TR(paramses["pi1.old.single"],cspan=c(1),se=TRUE,dec=3)%:%
    TR(c(newguess["pi1.old.married"]),cspan=c(1),dec=3)%:%TR(paramses["pi1.old.married"],cspan=c(1),se=TRUE,dec=3)%:%
    TR("---",cspan=2)+
    TR(c("$\\pi_2^{old}$"))%:%TR(newguess["pi2.old.single"],cspan=c(1),dec=3)%:%TR(paramses["pi2.old.single"],cspan=c(1),se=TRUE,dec=3)%:%
    TR(c(newguess["pi2.old.married"]),cspan=c(1),dec=3)%:%TR(paramses["pi2.old.married"],cspan=c(1),se=TRUE,dec=3)%:%
    TR("---",cspan=2)
  
  TS(printparams,file=paste0("paramvals_all_",modeltype),
     output_path=path,pretty_rules=T,
     header=c("r",rep("c",6)))
  
  
}
newguess




printwages<-TR(c("","Worker","Spouse"),cspan=c(1,2,2))+
  TR(c("","Coef","St. Error","Coef","St. Error"))+
  midrulep(list(c(2,3),c(4,5))) +
  TR("Panel A: Wage Coefficients",cspan=c(3), surround = "\\textit{%s}")+
  TR(c("Age"))%:%TR(c(newguess[c("age")]),cspan=c(1),dec=3)%:%TR(c(paramses[c("age")]),cspan=c(1),se=TRUE,dec=3)%:%
  TR(newguess["spage"],cspan=1,dec=3)%:%TR(paramses["spage"],cspan=1,se=TRUE,dec=3)+
  TR(c("Age$^2$"))%:%TR(c(newguess["age2"]),cspan=c(1),dec=3)%:%TR(paramses["age2"],cspan=c(1),se=TRUE,dec=3)%:%
  TR(newguess["spage2"],cspan=1,dec=3)%:%TR(paramses["spage2"],cspan=1,se=TRUE,dec=3)+
  TR("$H^{mod}$")%:%TR(newguess[c("mod")],cspan=c(1),dec=3)%:%TR(c(paramses[c("mod")]),cspan=c(1),dec=3,se=TRUE)%:%
  TR(c(newguess["spmod"]),cspan=c(1),dec=3)%:%TR(paramses["spmod"],cspan=c(1),se=TRUE,dec=3)+
  TR("$H^{sev}$")%:%TR(newguess[c("sev")],cspan=c(1),dec=3)%:%TR(c(paramses[c("sev")]),cspan=c(1),dec=3,se=TRUE)%:%
  TR("---",cspan=2)+
  vspace(5)+
  TR("Panel B: Type Effects",cspan=c(3), surround = "\\textit{%s}")+
  TR("$f_1$")%:%TR(newguess[c("Type1")],cspan=c(1),dec=3)%:%TR(c(paramses[c("Type1")]),cspan=c(1),dec=3,se=TRUE)%:%
  TR(c(newguess["spType1"]),cspan=c(1),dec=3)%:%TR(paramses["spType1"],cspan=c(1),se=TRUE,dec=3)+
  TR("$f_2$")%:%TR(newguess[c("Type2")],cspan=c(1),dec=3)%:%TR(c(paramses[c("Type2")]),cspan=c(1),dec=3,se=TRUE)%:%
  TR(c(newguess["spType2"]),cspan=c(1),dec=3)%:%TR(paramses["spType2"],cspan=c(1),se=TRUE,dec=3)+
  TR("$f_3$")%:%TR(newguess[c("Type3")],cspan=c(1),dec=3)%:%TR(c(paramses[c("Type3")]),cspan=c(1),dec=3,se=TRUE)%:%
  TR(c(newguess["spType3"]),cspan=c(1),dec=3)%:%TR(paramses["spType3"],cspan=c(1),se=TRUE,dec=3)+
  vspace(5)+
  TR("Panel C: Variance Parameters",cspan=c(3), surround = "\\textit{%s}")+
  TR("$\\sigma_\\epsilon$")%:%TR(c(newguess["noisevar"]^0.5),dec=4)%:%TR(c(paramses["noisevar"]^0.5), dec=4,se=TRUE)%:%
  TR(c(newguess["spnoisevar"]^0.5),dec=4)%:%TR(c(paramses["spnoisevar"]^0.5), dec=4,se=TRUE)+
  TR("$\\sigma_\\xi$")%:%TR(c(newguess["wagevar"]^0.5),dec=4)%:%TR(c(paramses["wagevar"]^0.5), dec=4,se=TRUE)%:%
  TR(c(newguess["spwagevar"]^0.5),dec=4)%:%TR(c(paramses["spwagevar"]^0.5), dec=4,se=TRUE)+
  TR("$\\rho_\\xi$")%:%TR(c(newguess["wagecorr"]),cspan=1,dec=4)%:%TR(c(paramses["wagecorr"]),cspan=1,dec=4,se=TRUE)%:%
  TR("---",cspan=2)
  


TS(printwages,file=paste0("paramvals_wages_",modeltype),
   output_path=path,pretty_rules=T,
   header=c("r",rep("c",4)))

# paramvals_A<-readRDS("../harmonized_code/modeloutput/paramvals_A.RDS")
# TS(paramvals_A,file="paramvals_A",
#    output_path=path,pretty_rules=T,
#    header=c("r",rep("c",6)))
# paramvals_B<-readRDS("../harmonized_code/modeloutput/paramvals_B.RDS")
# TS(paramvals_B,file="paramvals_B",
#    output_path=path,pretty_rules=T,
#    header=c("r",rep("c",6)))
# paramvals_C<-readRDS("../harmonized_code/modeloutput/paramvals_C.RDS")
# TS(paramvals_C,file="paramvals_C",
#    output_path=path,pretty_rules=T,
#    header=c("r",rep("c",6)))

wagecons<-readRDS(paste0("./modeloutput/momenttable_wagecons_",modeltype,".RDS"))
TS(wagecons,file=paste0("momenttable_wagecons_",modeltype),
   output_path=path,pretty_rules=T,
   header=c("r",rep("c",6)))

workDI<-readRDS(paste0("./modeloutput/momenttable_workDI_",modeltype,".RDS"))
TS(workDI,file=paste0("momenttable_workDI_",modeltype),
   output_path=path,pretty_rules=T,
   header=c("r",rep("c",6)))

jointwork<-readRDS(paste0("./modeloutput/jointwork_",modeltype,".RDS"))
TS(jointwork,file=paste0("jointwork_",modeltype),
   output_path=path,pretty_rules=T,
   header=c("r",rep("c",6)))

stockDIfit<-readRDS(paste0("./modeloutput/stockDIfit_",modeltype,".RDS"))
TS(stockDIfit,file=paste0("stockDIfit_",modeltype),
   output_path=path,pretty_rules=T,
   header=c("r",rep("c",4)))


wagefit<-readRDS(paste0("./modeloutput/wagefit_",modeltype,".RDS"))
TS(wagefit,file=paste0("wagefit_",modeltype),
   output_path=path,pretty_rules=T,
   header=c("r",rep("c",4)))

spwagefit<-readRDS(paste0("./modeloutput/spwagefit_",modeltype,".RDS"))
TS(spwagefit,file=paste0("spwagefit_",modeltype),
   output_path=path,pretty_rules=T,
   header=c("r",rep("c",4)))

workfit<-readRDS(paste0("./modeloutput/workfit_",modeltype,".RDS"))
TS(workfit,file=paste0("workfit_",modeltype),
   output_path=path,pretty_rules=T,
   header=c("r",rep("c",4)))

workfit<-readRDS(paste0("./modeloutput/DIfit_",modeltype,".RDS"))
TS(workfit,file=paste0("DIfit_",modeltype),
   output_path=path,pretty_rules=T,
   header=c("r",rep("c",4)))

consfit<-readRDS(paste0("./modeloutput/consfit_",modeltype,".RDS"))
TS(consfit,file=paste0("consfit_",modeltype),
   output_path=path,pretty_rules=T,
   header=c("r",rep("c",4)))

eventfit<-readRDS(paste0("./modeloutput/eventstudyheadsev_sim_healthresults_",modeltype,".RDS"))
TS(eventfit,file=paste0("eventstudyheadsev_sim_",modeltype),
   output_path=path,pretty_rules=T,
   header=c("r",rep("c",8)))

kfit<-readRDS(paste0("./modeloutput/kmeans/kmeansim_postperformance",modeltype,".RDS"))
#kfit<-cbind(rep(0,4),kfit)
colnames(kfit)[1]<-"1"
kfit<-kfit[!is.na(rownames(kfit)),] #missing type_stationary means out of sample
kfit<-kfit/sum(kfit)
kfit[,1]<-kfit[,1]+kfit[,is.na(colnames(kfit))] #missing type_propdpost  means should be assigned type 1
kfit<-kfit[,!is.na(colnames(kfit))]

kfit_tbl<-TR(c("Household Type"),cspan=c(4))+
  TR(c("K-Means","Model Probability"),cspan=c(1,3))+
  midrulep(list(c(1,1),c(2,4)))+
  TR(c("","Low","Mid","High"))+
  TR("Low")%:%TR(kfit[1,],dec=2)+
  TR("Mid")%:%TR(kfit[2,],dec=2)+
  TR("High")%:%TR(kfit[3,],dec=2)

TS(kfit_tbl,file="kmeans_postperformance",
   output_path=path,pretty_rules=T,
   header=c("r",rep("c",3)))


elasticityfit_all<-readRDS(paste0("./modeloutput/elasticityfit_all_",modeltype,".RDS"))


TS(elasticityfit_all,file=paste0("elasticityfit_all_",modeltype),
   output_path=path, pretty_rules = T,
   header=c('r',rep('c',8)))


output<-readRDS("./modeloutput/wtp_decompose.RDS")  

uflow2<-uflow/exp(newguess["theta_sev"]/2)
#CV for mod dis
uflow-uflow2

meanC<-simdata[age<62,mean(Ceq)]
meanC2<-((meanC^(-params$gamma))*(exp(newguess["eta"])/exp(newguess["theta_sev"]))^(1-params$gamma))^(-1/params$gamma)
#CV:
meanC2-meanC


modeltype <- "simplified"
cfact_vars<-list()
cfacts<-c("ss","fs","stringency","singstringency","marstringency","earlydi","meanstest","meanstest20k","meanstest40k","meanstest90k","ss_meanstest","fs_meanstest","stringency_meanstest","singstringency_meanstest","marstringency_meanstest","earlydi_meanstest")
vcovnames<-c("wtp","wtpcash","constax","constaxcash",
             "gain", "mvpf", 
             "work","work.m.1","work.m.2","work.m.3","work.s.1","work.s.2","work.s.3",
             "DI","DI.m.1","DI.m.2","DI.m.3","DI.s.1","DI.s.2","DI.s.3")

for(cfact in cfacts){
  temp<-tryCatch(unique(as.data.table(read.csv(file=paste0("./modeloutput/bootstrap/bootstrap_",cfact,"_",modeltype,".csv"),header=FALSE))),error=function(x) NA)
#  if(cfact%in%c("ss_meanstest","fs_meanstest","stringency_meanstest","singstringency_meanstest","marstringency_meanstest","earlydi_meanstest")){
#    temp2<-tryCatch(read.csv(file=paste0("./modeloutput/gradient_","meanstest","_",modeltype,".csv"),header=FALSE),error=function(x) NA)
#    temp<-temp-temp2
#  }
  if(any(is.na(temp))){
    cfact_vars[[length(cfact_vars)+1]]<-as.data.table(matrix(rep(NA,2*length(cfact_vars[[1]])),nrow=2))
  }
  temp<-temp[1:200,] #just so I can say we used exactly 200 replications.
  cfact_vars[[length(cfact_vars)+1]]<-temp
  rm(temp)
  colnames(cfact_vars[[length(cfact_vars)]])<-c(vcovnames,"particle")
  cfact_vars[[length(cfact_vars)]][,mvpf:=-mvpf]
  cfact_vars[[length(cfact_vars)]][,constaxcash:=-constaxcash]
  cfact_vars[[length(cfact_vars)]][,gain:=wtpcash-constaxcash]
  cfact_vars[[length(cfact_vars)]][constaxcash<0&wtpcash>=0,mvpf:=Inf]
}
names(cfact_vars)<-c("SSgenerosity","FSgenerosity","stringency","singstringency","marstringency","EarlyDI","MeansTest_1","MeansTest_20000","MeansTest_40000","MeansTest_90000","SSgenerositymeanstest_1","FSgenerositymeanstest_1","stringencymeanstest_1","singstringencymeanstest_1","marstringencymeanstest_1","EarlyDImeanstest_1")
rm(cfacts)



trimfun<-function(xdt) {
  xdt<-xdt[,!grepl("\\.x",names(xdt)),with=FALSE]
  names(xdt) <- gsub("\\.y","",names(xdt))
  return(xdt)
}
output<-lapply(output,trimfun)



setwd("./modeloutput/")
output$summary_origcompare[,margin:=margin/100]
output$summary_typecompare[,margin:=margin/100]
output$summary_marcompare[,margin:=margin/100]
output$summary_martypecompare[,margin:=margin/100]


###########################################################

#defining proportional reforms
for(reform in c("wtp_SSgenerosity","wtp_FSgenerosity","wtp_stringency","wtp_marstringency","wtp_EarlyDI",
                "wtp_SSgenerositymeanstest.1","wtp_stringencymeanstest.1","wtp_marstringencymeanstest.1","wtp_EarlyDImeanstest.1"
                #"wtp_SSgenerositymeanstest.20000","wtp_stringencymeanstest.20000","wtp_marstringencymeanstest.20000","wtp_EarlyDImeanstest.20000",
                #"wtp_SSgenerositymeanstest.40000","wtp_stringencymeanstest.40000","wtp_marstringencymeanstest.40000","wtp_EarlyDImeanstest.40000",
                #"wtp_SSgenerositymeanstest.90000","wtp_stringencymeanstest.90000","wtp_marstringencymeanstest.90000","wtp_EarlyDImeanstest.90000",
                #"wtp_SSgenerositymeanstestnoSSDI","wtp_stringencymeanstestnoSSDI","wtp_marstringencymeanstestnoSSDI","wtp_EarlyDImeanstestnoSSDI"
                )){
  #,"wtp_WorkApply" ,"wtp_WorkApplymeanstest" ,"wtp_WorkApplymeanstestnoSSDI"
  #"wtp_WorkDI"
  for(cfact in unique(output$summary_origcompare$counterfactual)){
    for(pre in c("")){
      output$summary_origcompare[counterfactual==cfact & dataset == "Counterfactual",eval(paste0("exanteprop",pre,"_",reform)):=get(paste0("exante",pre,"_",reform))/(
        output$summary_origcompare[counterfactual==cfact & dataset == "Baseline",get(paste0("exante",pre,"_",reform))]
      ) - 1]
      output$summary_origcompare[counterfactual==cfact & dataset == "Counterfactual"& is.na(get(paste0("exanteprop",pre,"_",reform))),eval(paste0("exanteprop",pre,"_",reform)):=0]
    }
    
    #proportional WTP for shawp components is relative to WTP for the total reform
    for(pre in c("all"#,"GV","RB","RW","IB","IM"
    )){
      if(pre=="all") pre2<-""
      else pre2<-pre
      output$summary_origcompare[counterfactual==cfact & dataset == "Baseline",eval(paste0("exanteprop",pre,"_",reform)):=get(paste0("exante_",pre2,reform))/(
        output$summary_origcompare[counterfactual==cfact & dataset == "Baseline",get(paste0("exante","_",reform))]
      )]
      output$summary_origcompare[counterfactual==cfact & dataset == "Baseline"& is.na(get(paste0("exanteprop",pre,"_",reform))),eval(paste0("exanteprop",pre,"_",reform)):=0]
    }
    
  }
}

cfacts<-unique(output$summary_origcompare$counterfactual)
reformname<-c("DI Generosity", "Food Stamp Generosity","Allowance Probability")
####################################################
#######SUMMARY FIGURES AND DECOMPOSITIONS OF EX-ANTE WTP
###################################################


print_mvpfs<-function(output,filename,cfact,baseval="Baseline",se){
  summary_mvpf<- output$summary_origcompare[counterfactual==cfact&margin==1.1,.(
    mvpf.SSgenerosity=mean(exante_wtp_SSgenerosity/constax_wtp_SSgenerosity,na.rm=TRUE),
    mvpf.SSgenerositymeanstest.1=mean(exante_wtp_SSgenerositymeanstest.1/constax_wtp_SSgenerositymeanstest.1,na.rm=TRUE),
    #mvpf.SSgenerositymeanstest.20000=mean(exante_wtp_SSgenerositymeanstest.20000/constax_wtp_SSgenerositymeanstest.20000,na.rm=TRUE),
    #      mvpf.SSgenerositymeanstest.40000=mean(exante_wtp_SSgenerositymeanstest.40000/constax_wtp_SSgenerositymeanstest.40000,na.rm=TRUE),
    #      mvpf.SSgenerositymeanstest.90000=mean(exante_wtp_SSgenerositymeanstest.90000/constax_wtp_SSgenerositymeanstest.90000,na.rm=TRUE),
    #mvpf.SSgenerositymeanstestnoSSDI=mean(exante_wtp_SSgenerositymeanstestnoSSDI/constax_wtp_SSgenerositymeanstestnoSSDI,na.rm=TRUE),
    mvpf.stringency=mean(exante_wtp_stringency/constax_wtp_stringency,na.rm=TRUE),
    mvpf.stringencymeanstest.1=mean(exante_wtp_stringencymeanstest.1/constax_wtp_stringencymeanstest.1,na.rm=TRUE),
    #mvpf.stringencymeanstest.20000=mean(exante_wtp_stringencymeanstest.20000/constax_wtp_stringencymeanstest.20000,na.rm=TRUE),
    #     mvpf.stringencymeanstest.40000=mean(exante_wtp_stringencymeanstest.40000/constax_wtp_stringencymeanstest.40000,na.rm=TRUE),
    #      mvpf.stringencymeanstest.90000=mean(exante_wtp_stringencymeanstest.90000/constax_wtp_stringencymeanstest.90000,na.rm=TRUE),
   # mvpf.stringencymeanstestnoSSDI=mean(exante_wtp_stringencymeanstestnoSSDI/constax_wtp_stringencymeanstestnoSSDI,na.rm=TRUE),
    mvpf.marstringency=mean(exante_wtp_marstringency/constax_wtp_marstringency,na.rm=TRUE),
    mvpf.marstringencymeanstest.1=mean(exante_wtp_marstringencymeanstest.1/constax_wtp_marstringencymeanstest.1,na.rm=TRUE),
    #mvpf.marstringencymeanstest.20000=mean(exante_wtp_marstringencymeanstest.20000/constax_wtp_marstringencymeanstest.20000,na.rm=TRUE),
    #     mvpf.marstringencymeanstest.40000=mean(exante_wtp_marstringencymeanstest.40000/constax_wtp_marstringencymeanstest.40000,na.rm=TRUE),
    #      mvpf.marstringencymeanstest.90000=mean(exante_wtp_marstringencymeanstest.90000/constax_wtp_marstringencymeanstest.90000,na.rm=TRUE),
    #mvpf.marstringencymeanstestnoSSDI=mean(exante_wtp_singstringencymeanstestnoSSDI/constax_wtp_singstringencymeanstestnoSSDI,na.rm=TRUE),
    mvpf.singstringency=mean(exante_wtp_singstringency/constax_wtp_singstringency,na.rm=TRUE),
    mvpf.singstringencymeanstest.1=mean(exante_wtp_singstringencymeanstest.1/constax_wtp_singstringencymeanstest.1,na.rm=TRUE),
    #mvpf.singstringencymeanstest.20000=mean(exante_wtp_singstringencymeanstest.20000/constax_wtp_singstringencymeanstest.20000,na.rm=TRUE),
    #     mvpf.singstringencymeanstest.40000=mean(exante_wtp_marstringencymeanstest.40000/constax_wtp_marstringencymeanstest.40000,na.rm=TRUE),
    #      mvpf.singstringencymeanstest.90000=mean(exante_wtp_marstringencymeanstest.90000/constax_wtp_marstringencymeanstest.90000,na.rm=TRUE),
    #mvpf.singstringencymeanstestnoSSDI=mean(exante_wtp_singstringencymeanstestnoSSDI/constax_wtp_singstringencymeanstestnoSSDI,na.rm=TRUE),
    mvpf.FSgenerosity=mean(exante_wtp_FSgenerosity/constax_wtp_FSgenerosity,na.rm=TRUE),
    mvpf.FSgenerositymeanstest.1=mean(exante_wtp_FSgenerositymeanstest.1/constax_wtp_FSgenerositymeanstest.1,na.rm=TRUE),
    #mvpf.FSgenerositymeanstest.20000=mean(exante_wtp_FSgenerositymeanstest.20000/constax_wtp_FSgenerositymeanstest.20000,na.rm=TRUE),
    #      mvpf.FSgenerositymeanstest.40000=mean(exante_wtp_FSgenerositymeanstest.40000/constax_wtp_FSgenerositymeanstest.40000,na.rm=TRUE),
    #      mvpf.FSgenerositymeanstest.90000=mean(exante_wtp_FSgenerositymeanstest.90000/constax_wtp_FSgenerositymeanstest.90000,na.rm=TRUE),
    #mvpf.FSgenerositymeanstestnoSSDI=mean(exante_wtp_FSgenerositymeanstestnoSSDI/constax_wtp_FSgenerositymeanstestnoSSDI,na.rm=TRUE),
    #      mvpf.WorkApply=mean(exante_wtp_WorkApply/constax_wtp_WorkApply,na.rm=TRUE),
    #      mvpf.WorkApplymeanstest=mean(exante_wtp_WorkApplymeanstest/constax_wtp_WorkApplymeanstest,na.rm=TRUE),
    #      mvpf.WorkApplymeanstestnoSSDI=mean(exante_wtp_WorkApplymeanstestnoSSDI/constax_wtp_WorkApplymeanstestnoSSDI,na.rm=TRUE),
    #      mvpf.WorkDI=mean(exante_wtp_WorkDI/constax_wtp_WorkDI,na.rm=TRUE),
    #      mvpf.WorkDImeanstest=mean(exante_wtp_WorkDImeanstest/constax_wtp_WorkDImeanstest,na.rm=TRUE),
    #      mvpf.WorkDImeanstestnoSSDI=mean(exante_wtp_WorkDImeanstestnoSSDI/constax_wtp_WorkDImeanstestnoSSDI,na.rm=TRUE),
    mvpf.EarlyDI=mean(exante_wtp_EarlyDI/constax_wtp_EarlyDI,na.rm=TRUE),
    mvpf.EarlyDImeanstest.1=mean(exante_wtp_EarlyDImeanstest.1/constax_wtp_EarlyDImeanstest.1,na.rm=TRUE)
    #mvpf.EarlyDImeanstest.20000=mean(exante_wtp_EarlyDImeanstest.20000/constax_wtp_EarlyDImeanstest.20000,na.rm=TRUE),
    #      mvpf.EarlyDImeanstest.40000=mean(exante_wtp_EarlyDImeanstest.40000/constax_wtp_EarlyDImeanstest.40000,na.rm=TRUE),
    #      mvpf.EarlyDImeanstest.90000=mean(exante_wtp_EarlyDImeanstest.90000/constax_wtp_EarlyDImeanstest.90000,na.rm=TRUE),
    #mvpf.EarlyDImeanstestnoSSDI=mean(exante_wtp_EarlyDImeanstestnoSSDI/constax_wtp_EarlyDImeanstestnoSSDI,na.rm=TRUE)
  ),
  by=.(counterfactual,dataset)]
  
  
  summary_mvpf_bymartype<- output$summary_martypecompare[counterfactual==cfact&margin==1.1,.(
    mvpf.SSgenerosity=mean(exanteMarType_wtp_SSgenerosity/constax_wtp_SSgenerosity,na.rm=TRUE),
    mvpf.SSgenerositymeanstest.1=mean(exanteMarType_wtp_SSgenerositymeanstest.1/constax_wtp_SSgenerositymeanstest.1,na.rm=TRUE),
    #mvpf.SSgenerositymeanstest.20000=mean(exanteMarType_wtp_SSgenerositymeanstest.20000/constax_wtp_SSgenerositymeanstest.20000,na.rm=TRUE),
    #      mvpf.SSgenerositymeanstest.40000=mean(exanteMarType_wtp_SSgenerositymeanstest.40000/constax_wtp_SSgenerositymeanstest.40000,na.rm=TRUE),
    #      mvpf.SSgenerositymeanstest.90000=mean(exanteMarType_wtp_SSgenerositymeanstest.90000/constax_wtp_SSgenerositymeanstest.90000,na.rm=TRUE),
    # mvpf.SSgenerositymeanstestnoSSDI=mean(exanteMarType_wtp_SSgenerositymeanstestnoSSDI/constax_wtp_SSgenerositymeanstestnoSSDI,na.rm=TRUE),
    mvpf.stringency=mean(exanteMarType_wtp_stringency/constax_wtp_stringency,na.rm=TRUE),
    mvpf.stringencymeanstest.1=mean(exanteMarType_wtp_stringencymeanstest.1/constax_wtp_stringencymeanstest.1,na.rm=TRUE),
    #mvpf.stringencymeanstest.20000=mean(exanteMarType_wtp_stringencymeanstest.20000/constax_wtp_stringencymeanstest.20000,na.rm=TRUE),
    #     mvpf.stringencymeanstest.40000=mean(exanteMarType_wtp_stringencymeanstest.40000/constax_wtp_stringencymeanstest.40000,na.rm=TRUE),
    #      mvpf.stringencymeanstest.90000=mean(exanteMarType_wtp_stringencymeanstest.90000/constax_wtp_stringencymeanstest.90000,na.rm=TRUE),
    # mvpf.stringencymeanstestnoSSDI=mean(exanteMarType_wtp_stringencymeanstestnoSSDI/constax_wtp_stringencymeanstestnoSSDI,na.rm=TRUE),
    mvpf.marstringency=mean(exanteMarType_wtp_marstringency/constax_wtp_marstringency,na.rm=TRUE),
    mvpf.marstringencymeanstest.1=mean(exanteMarType_wtp_marstringencymeanstest.1/constax_wtp_marstringencymeanstest.1,na.rm=TRUE),
    #mvpf.marstringencymeanstest.20000=mean(exanteMarType_wtp_marstringencymeanstest.20000/constax_wtp_marstringencymeanstest.20000,na.rm=TRUE),
    #     mvpf.marstringencymeanstest.40000=mean(exanteMarType_wtp_marstringencymeanstest.40000/constax_wtp_marstringencymeanstest.40000,na.rm=TRUE),
    #      mvpf.marstringencymeanstest.90000=mean(exanteMarType_wtp_marstringencymeanstest.90000/constax_wtp_marstringencymeanstest.90000,na.rm=TRUE),
    # mvpf.marstringencymeanstestnoSSDI=mean(exanteMarType_wtp_marstringencymeanstestnoSSDI/constax_wtp_marstringencymeanstestnoSSDI,na.rm=TRUE),
    mvpf.singstringency=mean(exanteMarType_wtp_singstringency/constax_wtp_singstringency,na.rm=TRUE),
    mvpf.singstringencymeanstest.1=mean(exanteMarType_wtp_singstringencymeanstest.1/constax_wtp_singstringencymeanstest.1,na.rm=TRUE),
    #mvpf.singstringencymeanstest.20000=mean(exanteMarType_wtp_singstringencymeanstest.20000/constax_wtp_singstringencymeanstest.20000,na.rm=TRUE),
    #     mvpf.singstringencymeanstest.40000=mean(exanteMarType_wtp_stringencymeanstest.40000/constax_wtp_stringencymeanstest.40000,na.rm=TRUE),
    #      mvpf.singstringencymeanstest.90000=mean(exanteMarType_wtp_stringencymeanstest.90000/constax_wtp_stringencymeanstest.90000,na.rm=TRUE),
    # mvpf.singstringencymeanstestnoSSDI=mean(exanteMarType_wtp_singstringencymeanstestnoSSDI/constax_wtp_singstringencymeanstestnoSSDI,na.rm=TRUE),
    mvpf.FSgenerosity=mean(exanteMarType_wtp_FSgenerosity/constax_wtp_FSgenerosity,na.rm=TRUE),
    mvpf.FSgenerositymeanstest.1=mean(exanteMarType_wtp_FSgenerositymeanstest.1/constax_wtp_FSgenerositymeanstest.1,na.rm=TRUE),
    #mvpf.FSgenerositymeanstest.20000=mean(exanteMarType_wtp_FSgenerositymeanstest.20000/constax_wtp_FSgenerositymeanstest.20000,na.rm=TRUE),
    #     mvpf.FSgenerositymeanstest.40000=mean(exanteMarType_wtp_FSgenerositymeanstest.40000/constax_wtp_FSgenerositymeanstest.40000,na.rm=TRUE),
    #      mvpf.FSgenerositymeanstest.90000=mean(exanteMarType_wtp_FSgenerositymeanstest.90000/constax_wtp_FSgenerositymeanstest.90000,na.rm=TRUE),
    # mvpf.FSgenerositymeanstestnoSSDI=mean(exanteMarType_wtp_FSgenerositymeanstestnoSSDI/constax_wtp_FSgenerositymeanstestnoSSDI,na.rm=TRUE),
    #      mvpf.WorkApply=mean(exanteMarType_wtp_WorkApply/constax_wtp_WorkApply,na.rm=TRUE),
    #      mvpf.WorkApplymeanstest=mean(exanteMarType_wtp_WorkApplymeanstest/constax_wtp_WorkApplymeanstest,na.rm=TRUE),
    #      mvpf.WorkApplymeanstestnoSSDI=mean(exanteMarType_wtp_WorkApplymeanstestnoSSDI/constax_wtp_WorkApplymeanstestnoSSDI,na.rm=TRUE),
    #      mvpf.WorkDI=mean(exanteMarType_wtp_WorkDI/constax_wtp_WorkDI,na.rm=TRUE),
    #      mvpf.WorkDImeanstest=mean(exanteMarType_wtp_WorkDImeanstest/constax_wtp_WorkDImeanstest,na.rm=TRUE),
    #      mvpf.WorkDImeanstestnoSSDI=mean(exanteMarType_wtp_WorkDImeanstestnoSSDI/constax_wtp_WorkDImeanstestnoSSDI,na.rm=TRUE),
    mvpf.EarlyDI=mean(exanteMarType_wtp_EarlyDI/constax_wtp_EarlyDI,na.rm=TRUE),
    mvpf.EarlyDImeanstest.1=mean(exanteMarType_wtp_EarlyDImeanstest.1/constax_wtp_EarlyDImeanstest.1,na.rm=TRUE)
    #mvpf.EarlyDImeanstest.20000=mean(exanteMarType_wtp_EarlyDImeanstest.20000/constax_wtp_EarlyDImeanstest.20000,na.rm=TRUE),
    #      mvpf.EarlyDImeanstest.40000=mean(exanteMarType_wtp_EarlyDImeanstest.40000/constax_wtp_EarlyDImeanstest.40000,na.rm=TRUE),
    #      mvpf.EarlyDImeanstest.90000=mean(exanteMarType_wtp_EarlyDImeanstest.90000/constax_wtp_EarlyDImeanstest.90000,na.rm=TRUE),
    #mvpf.EarlyDImeanstestnoSSDI=mean(exanteMarType_wtp_EarlyDImeanstestnoSSDI/constax_wtp_EarlyDImeanstestnoSSDI,na.rm=TRUE),
  ),
  by=.(counterfactual,dataset,Type,single)]
  
  summary_wtp<-output$summary_origcompare[counterfactual==cfact,]
  names(summary_wtp)<-gsub("wtp_","wtp.",names(summary_wtp))
  names(summary_wtp)<-gsub("meanstest\\.","meanstest_",names(summary_wtp))
  names(summary_wtp)<-gsub("MeansTest\\.","MeansTest_",names(summary_wtp))
  summary_wtp<-summary_wtp[,names(summary_wtp)%in%c("margin","dataset","counterfactual")|
                             grepl("emp_",names(summary_wtp)) | 
                             grepl("rolls_",names(summary_wtp)) | 
                             grepl("apply_",names(summary_wtp)) | 
                             grepl("exante_wtp",names(summary_wtp)) |
                             grepl("constaxcash_wtp",names(summary_wtp))|
                             grepl("exantecash_wtp",names(summary_wtp))|
                             grepl("selection_wtp",names(summary_wtp))|
                             grepl("moralhazard_wtp",names(summary_wtp))|
                             str_starts(names(summary_wtp),"constax_wtp"),with=FALSE]
  summary_wtp<-reshape(summary_wtp,
                       varying=names(summary_wtp)[!names(summary_wtp)%in%c("margin","dataset","counterfactual")],
                       direction="long",sep=".")
  setnames(summary_wtp,"time","reform")
  
  summary_wtp<-summary_wtp[((margin==1.1)),]
  
  summary_wtpbytype<-output$summary_typecompare[counterfactual==cfact,]
  names(summary_wtpbytype)<-gsub("wtp_","wtp.",names(summary_wtpbytype))
  names(summary_wtpbytype)<-gsub("meanstest\\.","meanstest_",names(summary_wtpbytype))
  names(summary_wtpbytype)<-gsub("MeansTest\\.","MeansTest_",names(summary_wtpbytype))
  summary_wtpbytype<-summary_wtpbytype[,names(summary_wtpbytype)%in%c("margin","dataset","counterfactual","Type")|
                             grepl("empType_",names(summary_wtpbytype)) | 
                             grepl("rollsType_",names(summary_wtpbytype)) | 
                             grepl("applyType_",names(summary_wtpbytype)) | 
                             grepl("exanteType_wtp",names(summary_wtpbytype)) |
                               grepl("constaxcash_wtp",names(summary_wtpbytype))|
                               grepl("constaxcashType_wtp",names(summary_wtpbytype))|
                               grepl("exantecash_wtp",names(summary_wtpbytype))|
                               grepl("exantecashType_wtp",names(summary_wtpbytype))|
                               str_starts(names(summary_wtpbytype),"constax_wtp"),with=FALSE]
  summary_wtpbytype<-reshape(summary_wtpbytype,
                       varying=names(summary_wtpbytype)[!names(summary_wtpbytype)%in%c("margin","dataset","counterfactual","Type")],
                       direction="long",sep=".")
  setnames(summary_wtpbytype,"time","reform")
  
  summary_wtpbytype<-summary_wtpbytype[((margin==1.1)),]
  
  summary_wtpbymar<-output$summary_marcompare[counterfactual==cfact,]
  names(summary_wtpbymar)<-gsub("wtp_","wtp.",names(summary_wtpbymar))
  names(summary_wtpbymar)<-gsub("meanstest\\.","meanstest_",names(summary_wtpbymar))
  names(summary_wtpbymar)<-gsub("MeansTest\\.","MeansTest_",names(summary_wtpbymar))
  summary_wtpbymar<-summary_wtpbymar[,names(summary_wtpbymar)%in%c("margin","dataset","counterfactual","single")|
                                         grepl("rollsMar",names(summary_wtpbymar)),with=FALSE]
  summary_wtpbymar<-reshape(summary_wtpbymar,
                             varying=names(summary_wtpbymar)[!names(summary_wtpbymar)%in%c("margin","dataset","counterfactual","single")],
                             direction="long",sep=".")
  setnames(summary_wtpbymar,"time","reform")
  
  summary_wtpbymar<-summary_wtpbymar[((margin==1.1)),]
  
  
  summary_wtpbymartype<-output$summary_martypecompare[counterfactual==cfact,]
  names(summary_wtpbymartype)<-gsub("wtp_","wtp.",names(summary_wtpbymartype))
  names(summary_wtpbymartype)<-gsub("meanstest\\.","meanstest_",names(summary_wtpbymartype))
  names(summary_wtpbymartype)<-gsub("MeansTest\\.","MeansTest_",names(summary_wtpbymartype))
  summary_wtpbymartype<-summary_wtpbymartype[,names(summary_wtpbymartype)%in%c("margin","dataset","counterfactual","Type","single")|
                                               grepl("empMarType_",names(summary_wtpbymartype)) | 
                                               grepl("rollsMarType",names(summary_wtpbymartype)) | 
                                               grepl("applyMarType_",names(summary_wtpbymartype)) | 
                                               grepl("everapplyMarType_",names(summary_wtpbymartype)) | 
                                               grepl("everapplyMar",names(summary_wtpbymartype)) | 
                                               grepl("exanteMarType_wtp",names(summary_wtpbymartype))|
                                               grepl("constaxcash_wtp",names(summary_wtpbymartype))|
                                               grepl("exantecash_wtp",names(summary_wtpbymartype))|
                                               grepl("selection",names(summary_wtpbymartype))|
                                               str_starts(names(summary_wtpbymartype),"constax_wtp"),with=FALSE]
  summary_wtpbymartype<-reshape(summary_wtpbymartype,
                                varying=names(summary_wtpbymartype)[!names(summary_wtpbymartype)%in%c("Type","single","margin","dataset","counterfactual")],
                                direction="long",sep=".")
  setnames(summary_wtpbymartype,"time","reform")
  
  summary_wtpbymartype<-summary_wtpbymartype[((margin==1.1)),]

  
  for(data in c("Baseline","Counterfactual")){
    if(se==TRUE & data=="Baseline"){
    tab<-TR(c("","(1)","(2)","(3)","(4)","(5)","(6)"))+
      TR(c("","Fiscal ","Ex-Ante","Ex-Ante","Prob.","Total","Beneficiary"
           #, "ex-post (approximate)"
      ))
      tab<-tab+TR(c("","Costs (\\$)","WTP (\\$)","MVPF", "(3)$\\ge1$","Emp. (p.p.)","Rolls (p.p.)"
           #,"WTP (\\$)","MVPF"
      ))+
      midrulep(list(c(2,7)))
    } else{
      tab<-TR(c("","(1)","(2)","(3)","(4)","(5)"))+
        TR(c("","Fiscal ","Ex-Ante","Ex-Ante","Total","Beneficiary"
             #, "ex-post (approximate)"
        ))
      tab<-tab+TR(c("","Costs (\\$)","WTP (\\$)","MVPF", "Emp. (p.p.)","Rolls (p.p.)"
                    #,"WTP (\\$)","MVPF"
      ))+
        midrulep(list(c(2,6)))
    }
    tab<- tab + TR("Panel A: Effect of Baseline Reforms",cspan=c(6), surround = "\\textit{%s}")
    for(r in which(reformnames%in%c("EarlyDI","stringency","singstringency","SSgenerosity","FSgenerosity"))){
      tabrow<-TR(reformlabels[r])%:%
        TR(unlist(summary_wtp[reform==reformnames[r] & dataset==data,constaxcash_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames[r] & dataset==data,exantecash_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames[r] & dataset==data,exante_wtp/constax_wtp]),dec=2)
        if(se == TRUE & data=="Baseline") tabrow<-tabrow%:% TR(unlist(mean(cfact_vars[[which(names(cfact_vars)==reformnames[r])]][,"mvpf"]>=1)),dec=2)
        tabrow<-tabrow%:%TR(unlist(summary_wtp[reform==reformnames[r] & dataset==data,emp_wtp*100]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames[r] & dataset==data,rolls_wtp*100]),dec=2)
      tab<-tab+tabrow
      if(se ==TRUE & data=="Baseline") {
        tab<-tab+
          TR(c(""))%:%TR(unlist(
            apply(cfact_vars[[which(names(cfact_vars)==reformnames[r])]][,c("constaxcash","wtpcash","mvpf")],
                  MARGIN=2,FUN=function(x) paste0(c("[",paste0(round(quantile(x,probs=c(0.025,0.975)),2),collapse=", "),"]"),collapse="")
                  )))%:%TR(c(""))%:%
          TR(unlist(
            apply(cfact_vars[[which(names(cfact_vars)==reformnames[r])]][,c("work","DI")],
                  MARGIN=2,FUN=function(x) paste0(c("[",paste0(round(quantile(x,probs=c(0.025,0.975)),2),collapse=", "),"]"),collapse="")
            )))
      }
    }
    
    tab<- tab + vspace(5)  + TR("Panel B: Marginal Effect, With Strict Asset-Tested DI",cspan=c(6), surround = "\\textit{%s}")
    for(r in which(reformnames_meanstest%in%c("EarlyDImeanstest_1","stringencymeanstest_1","singstringencymeanstest_1","SSgenerositymeanstest_1","FSgenerositymeanstest_1"))){
      tabrow<-TR(reformlabels_meanstest[r])%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstest[r] & dataset==data,constaxcash_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstest[r] & dataset==data,exantecash_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstest[r] & dataset==data,exante_wtp/constax_wtp]),dec=2)
        if(se==TRUE & data=="Baseline") tabrow<-tabrow%:%TR(unlist(mean(cfact_vars[[which(names(cfact_vars)==reformnames_meanstest[r])]][,"mvpf"]>=1)),dec=2)
        tabrow<-tabrow%:%TR(unlist(summary_wtp[reform==reformnames_meanstest[r] & dataset==data,emp_wtp*100]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstest[r] & dataset==data,rolls_wtp*100]),dec=2)
        tab<-tab+tabrow
      if(se ==TRUE & data=="Baseline") {
        tab<-tab+
          TR(c(""))%:%TR(unlist(
            apply(cfact_vars[[which(names(cfact_vars)==reformnames_meanstest[r])]][,c("constaxcash","wtpcash","mvpf")],
                  MARGIN=2,FUN=function(x) paste0(c("[",paste0(round(quantile(x,probs=c(0.025,0.975)),2),collapse=", "),"]"),collapse="")
            )))%:%TR(c(""))%:%
          TR(unlist(
            apply(cfact_vars[[which(names(cfact_vars)==reformnames_meanstest[r])]][,c("work","DI")],
                  MARGIN=2,FUN=function(x) paste0(c("[",paste0(round(quantile(x,probs=c(0.025,0.975)),2),collapse=", "),"]"),collapse="")
            )))
      }
    }
    TS(tab,file=paste0("mvpf_abridged",data,filename),
       output_path="./",
       pretty_rules=T,
       header=c('r',rep('c',5+(se==TRUE  & data=="Baseline"))))
    
    
    tab<-TR(c("","(1)","(2)","(3)","(4)","(5)"))+
      TR(c("","Fiscal ","Ex-Ante","Ex-Ante","Total","Beneficiary"
           #, "ex-post (approximate)"
      ))+
      TR(c("","Costs (\\$)","WTP (\\$)","MVPF","Emp. (p.p.)","Rolls (p.p.)"
           #,"WTP (\\$)","MVPF"
      ))+
      midrulep(list(c(2,6)))
    tab<- tab + TR("Panel A: Effect of Baseline Reforms",cspan=c(6), surround = "\\textit{%s}")
    for(r in which(reformnames%in%c("EarlyDI","stringency","singstringency","SSgenerosity","FSgenerosity"))){
      tab<-tab+TR(reformlabels[r])%:%
        TR(unlist(summary_wtp[reform==reformnames[r] & dataset==data,constaxcash_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames[r] & dataset==data,exantecash_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames[r] & dataset==data,exante_wtp/constax_wtp]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames[r] & dataset==data,emp_wtp*100]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames[r] & dataset==data,rolls_wtp*100]),dec=2)
    }
    
    tab<- tab + vspace(5) +
      TR("Panel B: Effect of Means Tests",cspan=c(6), surround = "\\textit{%s}")
    for(r in 1:length(reformnames_meanstestonly)){
      tab<-tab+TR(reformlabels_meanstestonly[r])%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstestonly[r] & dataset==data,constaxcash_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstestonly[r] & dataset==data,exantecash_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstestonly[r] & dataset==data,exante_wtp/constax_wtp]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstestonly[r] & dataset==data,emp_wtp*100]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstestonly[r] & dataset==data,rolls_wtp*100]),dec=2)
    }
    tab<- tab + vspace(5)  + TR("Panel C: Marginal Effect, With Strict Asset-Tested DI",cspan=c(6), surround = "\\textit{%s}")
    for(r in which(reformnames_meanstest%in%c("EarlyDImeanstest_1","stringencymeanstest_1","singstringencymeanstest_1","SSgenerositymeanstest_1","FSgenerositymeanstest_1"))){
      tab<-tab+TR(reformlabels_meanstest[r])%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstest[r] & dataset==data,constaxcash_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstest[r] & dataset==data,exantecash_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstest[r] & dataset==data,exante_wtp/constax_wtp]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstest[r] & dataset==data,emp_wtp*100]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstest[r] & dataset==data,rolls_wtp*100]),dec=2)

    }
    TS(tab,file=paste0("mvpf_combined",data,filename),
       output_path="./",
       pretty_rules=T,
       header=c('r',rep('c',5)))
    
#    if(data=="Baseline"){
    ##WTP BY TYPE:
    tab<-TR(c("","(1)","(2)","(3)","(4)","(5)","(6)","(7)","(8)"))+
      TR(c("","Overall","Low Type","Mid Type", "High Type"),cspan=c(1,2,2,2,2))+
      TR(c("","Costs (\\$)","WTP (\\$)","Costs (\\$)","WTP (\\$)","Costs (\\$)","WTP (\\$)","Costs (\\$)","WTP (\\$)"))+
      midrulep(list(c(2,3),c(4,5),c(6,7),c(8,9)))
    tab<- tab + TR("Panel A: Effect of Baseline Reforms",cspan=c(9), surround = "\\textit{%s}")
    for(r in which(reformnames%in%c("EarlyDI","stringency","singstringency","SSgenerosity","FSgenerosity"))){
      tab<-tab+TR(reformlabels[r])%:%
        TR(unlist(summary_wtp[reform==reformnames[r] & dataset==data,constaxcash_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames[r] & dataset==data,exantecash_wtp]),dec=2)%:%
        TR(unlist(ifelse(length(summary_wtpbytype[reform==reformnames[r] & dataset==data & Type==1,constaxcashType_wtp])==0,
                         NA,summary_wtpbytype[reform==reformnames[r] & dataset==data & Type==1,constaxcashType_wtp])),dec=2)%:%
        TR(unlist(ifelse(length(summary_wtpbytype[reform==reformnames[r] & dataset==data & Type==1,exantecashType_wtp])==0,
                         NA,summary_wtpbytype[reform==reformnames[r] & dataset==data & Type==1,exantecashType_wtp])),dec=2)%:%
        TR(unlist(summary_wtpbytype[reform==reformnames[r] & dataset==data & Type==2,constaxcashType_wtp]),dec=2)%:%
        TR(unlist(summary_wtpbytype[reform==reformnames[r] & dataset==data & Type==2,exantecashType_wtp]),dec=2)%:%
        TR(unlist(summary_wtpbytype[reform==reformnames[r] & dataset==data & Type==3,constaxcashType_wtp]),dec=2)%:%
        TR(unlist(summary_wtpbytype[reform==reformnames[r] & dataset==data & Type==3,exantecashType_wtp]),dec=2)
        # if(se ==TRUE) tab<-tab+
        #   TR(c(""))%:%TR(cfact_vars[[r]][c("constax","wtp","mvpf","work","DI")],dec=2,ses=TRUE)
    }
    tab<- tab + vspace(5) +
     TR("Panel B: Effect of Means Tests",cspan=c(9), surround = "\\textit{%s}")
    for(r in 1:length(reformnames_meanstestonly)){
      tab<-tab+TR(reformlabels_meanstestonly[r])%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstestonly[r] & dataset==data,constaxcash_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstestonly[r] & dataset==data,exantecash_wtp]),dec=2)%:%
        TR(unlist(ifelse(length(summary_wtpbytype[reform==reformnames_meanstestonly[r] & dataset==data & Type==1,constaxcashType_wtp])==0,
                         NA,summary_wtpbytype[reform==reformnames_meanstestonly[r] & dataset==data & Type==1,constaxcashType_wtp])),dec=2)%:%
        TR(unlist(ifelse(length(summary_wtpbytype[reform==reformnames_meanstestonly[r] & dataset==data & Type==1,exantecashType_wtp])==0,
                         NA,summary_wtpbytype[reform==reformnames_meanstestonly[r] & dataset==data & Type==1,exantecashType_wtp])),dec=2)%:%
        TR(unlist(summary_wtpbytype[reform==reformnames_meanstestonly[r] & dataset==data & Type==2,constaxcashType_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtpbytype[reform==reformnames_meanstestonly[r] & dataset==data & Type==2,exantecashType_wtp]),dec=2)%:%
        TR(unlist(summary_wtpbytype[reform==reformnames_meanstestonly[r] & dataset==data & Type==3,constaxcashType_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtpbytype[reform==reformnames_meanstestonly[r] & dataset==data & Type==3,exantecashType_wtp]),dec=2)
    }
    
    
    tab<- tab + vspace(5)  + TR("Panel C: Marginal Effect, With Strict Asset-Tested DI",cspan=c(9), surround = "\\textit{%s}")
    for(r in which(reformnames_meanstest%in%c("EarlyDImeanstest_1","stringencymeanstest_1","singstringencymeanstest_1","SSgenerositymeanstest_1","FSgenerositymeanstest_1"))){
      tab<-tab+TR(reformlabels_meanstest[r])%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstest[r] & dataset==data,constaxcash_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstest[r] & dataset==data,exantecash_wtp/1]),dec=2)%:%
        TR(unlist(ifelse(length(summary_wtpbytype[reform==reformnames_meanstest[r] & dataset==data & Type==1,constaxcashType_wtp])==0,
                         NA,summary_wtpbytype[reform==reformnames_meanstest[r] & dataset==data & Type==1,constaxcashType_wtp])),dec=2)%:%
        TR(unlist(ifelse(length(summary_wtpbytype[reform==reformnames_meanstest[r] & dataset==data & Type==1,exantecashType_wtp])==0,
                         NA,summary_wtpbytype[reform==reformnames_meanstest[r] & dataset==data & Type==1,exantecashType_wtp])),dec=2)%:%
        TR(unlist(summary_wtpbytype[reform==reformnames_meanstest[r] & dataset==data & Type==2,constaxcashType_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtpbytype[reform==reformnames_meanstest[r] & dataset==data & Type==2,exantecashType_wtp]),dec=2)%:%
        TR(unlist(summary_wtpbytype[reform==reformnames_meanstest[r] & dataset==data & Type==3,constaxcashType_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtpbytype[reform==reformnames_meanstest[r] & dataset==data & Type==3,exantecashType_wtp]),dec=2)
#      if(se ==TRUE) tab<-tab+
#         TR(c(""))%:%TR(cfact_vars[[r]][c("constax","wtp","mvpf","work","DI")],dec=2,ses=TRUE)
    }
    TS(tab,file=paste0("wtp_bytype",data,filename),
       output_path="./",
       pretty_rules=T,
       header=c('r',rep('c',8)))
#    }
    
  }
  
  
  
    
    

  
  
  
  summary_wtp[reform%in%reformnames,group:="Baseline Reforms"]
  summary_wtp[reform%in%reformnames_meanstestonly,group:="Means Test"]
  summary_wtp[reform%in%reformnames_meanstest,group:="After Means Test"]
  summary_wtp[,group:=factor(group,levels=c("Baseline Reforms","Means Test","After Means Test"))]
  allreforms<-c(reformlabels,c("No Savings","Savings \u2264 20,000","Savings \u2264 40,000","Savings \u2264 90,000"),reformlabels_meanstest)
  allreforms[c(5,14)]<-"Food Stamp Generosity"
  names(allreforms)<-c(reformnames,
                       reformnames_meanstestonly,
                       reformnames_meanstest)
  for(data in c("Baseline"#,"Counterfactual"
                )){
    
    cairo_pdf(paste0("selectiongraph_",data,filename,".pdf"))
    print(ggplot(data=summary_wtp[dataset==data&reform%in%c(reformnames_meanstestonly,
                                                            reformnames,
                                                            reformnames_meanstest),])+
            geom_col(aes(x=as.factor(reform),y=selection_wtp*100)) + 
            labs(x='',y="Percent Change in Average Value of DI \n to DI Applicants",fill="",linetype="",pattern="")
          #      +scale_linetype_manual(values=c("Give-away"="solid","Redist. Between"="solid","Redist. Within"="solid",
          #                                      "Insur. (Health)" = "dashed","Insur. (Singlehood)"="dashed",
          #                                      "Insur. (Productivity)"="dashed","Non-pecuniary Effects"="dashed"))
          +geom_hline(yintercept=0,linetype="dashed")
          +coord_flip(ylim=c(-15,30))
          +facet_grid(rows=vars(group), scales="free_y", space="free_y", switch="y"
                      #labeller=labeller(group=grouplabs)
          )+
            theme(strip.text.y=element_text(size=14, angle=270),
                  strip.background = element_blank(),
                  axis.title.x = element_text(size=14),
                  axis.title.y = element_text(size=14),
                  axis.text.x = element_text(size=12),
                  axis.text.y = element_text(size=12),
                  strip.placement="outside")+
            scale_x_discrete(labels=allreforms)
#          +guides(fill=guide_legend(nrow=3),linetype=guide_legend(nrow=3),pattern=guide_legend(nrow=3))    
          )
    dev.off()
    
    tab<-TR(c("","(1)","(2)","(3)","(4)","(5)","(6)"))+
      TR(c("Marital Status:","Single","Married"),cspan=c(1,3,3))+
      TR(c("Type:","Low","Mid","High","Low","Mid","High"))+
      midrulep(list(c(2,4),c(5,7)))
    tab<- tab + TR("Panel A: Effect of Baseline Reforms",cspan=c(7), surround = "\\textit{%s}")
    for(r in which(reformnames%in%c("EarlyDI","stringency","singstringency","SSgenerosity","FSgenerosity","Meanstest_1","Meanstest_20000","Meanstest_40000","Meanstest_90000"))){
      tab<-tab+TR(reformlabels[r])%:%
        TR(unlist(summary_wtpbymartype[reform==reformnames[r] & dataset==data & Type==1&single==1,selectionMarType_wtp*100]),dec=2)%:%
        TR(unlist(summary_wtpbymartype[reform==reformnames[r] & dataset==data & Type==2&single==1,selectionMarType_wtp*100]),dec=2)%:%
        TR(unlist(summary_wtpbymartype[reform==reformnames[r] & dataset==data & Type==3&single==1,selectionMarType_wtp*100]),dec=2)%:%
        TR(unlist(summary_wtpbymartype[reform==reformnames[r] & dataset==data & Type==1&single==0,selectionMarType_wtp*100]),dec=2)%:%
        TR(unlist(summary_wtpbymartype[reform==reformnames[r] & dataset==data & Type==2&single==0,selectionMarType_wtp*100]),dec=2)%:%
        TR(unlist(summary_wtpbymartype[reform==reformnames[r] & dataset==data & Type==3&single==0,selectionMarType_wtp*100]),dec=2)
    }
    tab<- tab + vspace(5) +
      TR("Panel B: Effect of Means Tests",cspan=c(6), surround = "\\textit{%s}")
    for(r in 1:length(reformnames_meanstestonly)){
      tab<-tab+TR(reformlabels_meanstestonly[r])%:%
        TR(unlist(summary_wtpbymartype[reform==reformnames_meanstestonly[r] & dataset==data& Type==1&single==0,selectionMarType_wtp]),dec=2)%:%
        TR(unlist(summary_wtpbymartype[reform==reformnames_meanstestonly[r] & dataset==data& Type==2&single==0,selectionMarType_wtp]),dec=2)%:%
        TR(unlist(summary_wtpbymartype[reform==reformnames_meanstestonly[r] & dataset==data& Type==3&single==0,selectionMarType_wtp]),dec=2)%:%
        TR(unlist(summary_wtpbymartype[reform==reformnames_meanstestonly[r] & dataset==data & Type==1&single==1,selectionMarType_wtp]),dec=2)%:%
        TR(unlist(summary_wtpbymartype[reform==reformnames_meanstestonly[r] & dataset==data & Type==2&single==1,selectionMarType_wtp]),dec=2)%:%
        TR(unlist(summary_wtpbymartype[reform==reformnames_meanstestonly[r] & dataset==data & Type==3&single==1,selectionMarType_wtp]),dec=2)
    }
    
    tab<- tab + vspace(5)  + TR("Panel C: Marginal Effect, With Strict Asset-Tested DI",cspan=c(7), surround = "\\textit{%s}")
    for(r in which(reformnames_meanstest%in%c("EarlyDImeanstest_1","stringencymeanstest_1","singstringencymeanstest_1","SSgenerositymeanstest_1","FSgenerositymeanstest_1"))){
      tab<-tab+TR(reformlabels_meanstest[r])%:%
        TR(unlist(summary_wtpbymartype[reform==reformnames_meanstest[r] & dataset==data & Type==1&single==1,selectionMarType_wtp*100]),dec=2)%:%
        TR(unlist(summary_wtpbymartype[reform==reformnames_meanstest[r] & dataset==data & Type==2&single==1,selectionMarType_wtp*100]),dec=2)%:%
        TR(unlist(summary_wtpbymartype[reform==reformnames_meanstest[r] & dataset==data & Type==3&single==1,selectionMarType_wtp*100]),dec=2)%:%
        TR(unlist(summary_wtpbymartype[reform==reformnames_meanstest[r] & dataset==data & Type==1&single==0,selectionMarType_wtp*100]),dec=2)%:%
        TR(unlist(summary_wtpbymartype[reform==reformnames_meanstest[r] & dataset==data & Type==2&single==0,selectionMarType_wtp*100]),dec=2)%:%
        TR(unlist(summary_wtpbymartype[reform==reformnames_meanstest[r] & dataset==data & Type==3&single==0,selectionMarType_wtp*100]),dec=2)     }
    TS(tab,file=paste0("DIselection_",data,filename),
       output_path="./",
       pretty_rules=T,
       header=c('r',rep('c',6)))
    
  }
  
  #Effects of asset tests:
  for(data in c("Baseline","Counterfactual")){
    if(se==TRUE & data=="Baseline"){
    tab<-TR(c("","(1)","(2)","(3)","(4)","(5)","(6)"))+
      TR(c("","Fiscal ","Ex-Ante","Gain","Prob.","Total","Beneficiary"
           #, "ex-post (approximate)"
      ))+
      TR(c("","Costs (\\$)","WTP (\\$)","((2) - (1))","(3)$\\ge$0","Emp. (p.p.)","Rolls (p.p.)"
           #,"WTP (\\$)","MVPF"
      ))+
      midrulep(list(c(2,7)))
    }
    else {
      tab<-TR(c("","(1)","(2)","(3)","(4)","(5)"))+
        TR(c("","Fiscal ","Ex-Ante","Gain","Total","Beneficiary"
             #, "ex-post (approximate)"
        ))+
        TR(c("","Costs (\\$)","WTP (\\$)","((2) - (1))","Emp. (p.p.)","Rolls (p.p.)"
             #,"WTP (\\$)","MVPF"
        ))+
        midrulep(list(c(2,6)))
    }
    
    for(r in 1:length(reformnames_meanstestonly)){
      tabrow<-TR(reformlabels_meanstestonly[r])%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstestonly[r] & dataset==data,constaxcash_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstestonly[r] & dataset==data,exantecash_wtp/1]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstestonly[r] & dataset==data,exantecash_wtp - constaxcash_wtp]),dec=2)
      if(se ==TRUE & data=="Baseline") {
        tabrow<-tabrow%:%TR(unlist(mean(cfact_vars[[which(names(cfact_vars)==reformnames_meanstestonly[r])]][,"gain"]>=0)),dec=2)
      }
        tabrow<-tabrow%:%TR(unlist(summary_wtp[reform==reformnames_meanstestonly[r] & dataset==data,emp_wtp*100]),dec=2)%:%
        TR(unlist(summary_wtp[reform==reformnames_meanstestonly[r] & dataset==data,rolls_wtp*100]),dec=2)
        tab<-tab+tabrow
      if(se ==TRUE & data=="Baseline") {
        tab<-tab+
          TR(c(""))%:%TR(unlist(
            apply(cfact_vars[[which(names(cfact_vars)==reformnames_meanstestonly[r])]][,c("constaxcash","wtpcash","gain")],
                  MARGIN=2,FUN=function(x) paste0(c("[",paste0(round(quantile(x,probs=c(0.025,0.975)),2),collapse=", "),"]"),collapse="")
            )))%:%TR(c(""))%:%TR(unlist(apply(cfact_vars[[which(names(cfact_vars)==reformnames_meanstestonly[r])]][,c("work","DI")],
                                    MARGIN=2,FUN=function(x) paste0(c("[",paste0(round(quantile(x,probs=c(0.025,0.975)),2),collapse=", "),"]"),collapse="")
            )))
      }
    }
    TS(tab,file=paste0("wtp_assettests",data,filename),
       output_path="./",
       pretty_rules=T,
       header=c('r',rep('c',5+(se==TRUE & data=="Baseline"))))
  }
  
  
  
  
  
}




pctdiff<-function(cfact, reformnames){
  100*(unlist(cfact)/unlist(output$summary_origcompare[counterfactual=="nospwork" & dataset=="Baseline",(paste0("exante_wtp_",reformnames[r],c("","meanstest"))),with=FALSE])-1)
}



setwd("./modeloutput")

for(cfact in c("notype1",#"statedep",
               "fixmar")){
  if(cfact != "notype1") se <- TRUE
  else se <- FALSE
print_mvpfs(output=output,filename=cfact,cfact=cfact,se=se)
}
